From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 1 Aug 2020 12:45:54 +0000 (+0200)
Subject: Use face_no in FE
X-Git-Tag: v9.3.0-rc1~1191^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5080b6916e227801a5b0e83983e26141846d387a;p=dealii.git

Use face_no in FE
---

diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h
index 9763abf10d..88b5787810 100644
--- a/include/deal.II/fe/fe.h
+++ b/include/deal.II/fe/fe.h
@@ -3129,9 +3129,12 @@ template <int dim, int spacedim>
 inline std::pair<unsigned int, unsigned int>
 FiniteElement<dim, spacedim>::face_system_to_component_index(
   const unsigned int index,
-  const unsigned int) const
+  const unsigned int face_no) const
 {
-  AssertIndexRange(index, face_system_to_component_table[0].size());
+  AssertIndexRange(
+    index,
+    face_system_to_component_table[this->n_unique_faces() == 1 ? 0 : face_no]
+      .size());
 
   // in debug mode, check whether the
   // function is primitive, since
@@ -3146,11 +3149,13 @@ FiniteElement<dim, spacedim>::face_system_to_component_index(
   //
   // in 1d, the face index is equal
   // to the cell index
-  Assert(is_primitive(this->face_to_cell_index(index, 0)),
+  Assert(is_primitive(this->face_to_cell_index(index, face_no)),
          (typename FiniteElement<dim, spacedim>::ExcShapeFunctionNotPrimitive(
            index)));
 
-  return face_system_to_component_table[0][index];
+  return face_system_to_component_table[this->n_unique_faces() == 1 ?
+                                          0 :
+                                          face_no][index];
 }
 
 
@@ -3170,10 +3175,14 @@ template <int dim, int spacedim>
 inline std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
 FiniteElement<dim, spacedim>::face_system_to_base_index(
   const unsigned int index,
-  const unsigned int) const
+  const unsigned int face_no) const
 {
-  AssertIndexRange(index, face_system_to_base_table[0].size());
-  return face_system_to_base_table[0][index];
+  AssertIndexRange(
+    index,
+    face_system_to_base_table[this->n_unique_faces() == 1 ? 0 : face_no]
+      .size());
+  return face_system_to_base_table[this->n_unique_faces() == 1 ? 0 : face_no]
+                                  [index];
 }
 
 
@@ -3301,7 +3310,7 @@ FiniteElement<dim, spacedim>::get_associated_geometry_primitive(
   // are enumerated on the reference cell
   if (cell_dof_index < this->get_first_line_index())
     return GeometryPrimitive::vertex;
-  else if (cell_dof_index < this->get_first_quad_index())
+  else if (cell_dof_index < this->get_first_quad_index(0))
     return GeometryPrimitive::line;
   else if (cell_dof_index < this->get_first_hex_index())
     return GeometryPrimitive::quad;
diff --git a/include/deal.II/fe/fe_poly_face.templates.h b/include/deal.II/fe/fe_poly_face.templates.h
index fb94219cc1..06bc3ed006 100644
--- a/include/deal.II/fe/fe_poly_face.templates.h
+++ b/include/deal.II/fe/fe_poly_face.templates.h
@@ -141,16 +141,16 @@ FE_PolyFace<PolynomialType, dim, spacedim>::fill_fe_face_values(
             case 3:
               {
                 // Fill data for quad shape functions
-                if (this->n_dofs_per_quad() != 0)
+                if (this->n_dofs_per_quad(face_no) != 0)
                   {
                     const unsigned int foffset =
-                      this->get_first_quad_index() +
-                      this->n_dofs_per_quad() * face_no;
-                    for (unsigned int k = 0; k < this->n_dofs_per_quad(); ++k)
+                      this->get_first_quad_index(face_no);
+                    for (unsigned int k = 0; k < this->n_dofs_per_quad(face_no);
+                         ++k)
                       output_data.shape_values(foffset + k, i) =
                         fe_data
-                          .shape_values[k + this->get_first_face_quad_index()]
-                                       [i];
+                          .shape_values[k + this->get_first_face_quad_index(
+                                              face_no)][i];
                   }
               }
               DEAL_II_FALLTHROUGH;
@@ -173,9 +173,10 @@ FE_PolyFace<PolynomialType, dim, spacedim>::fill_fe_face_values(
                                                                     line) *
                                 this->n_dofs_per_line() +
                               k,
-                            i) = fe_data.shape_values
-                                   [k + (line * this->n_dofs_per_line()) +
-                                    this->get_first_face_line_index()][i];
+                            i) =
+                            fe_data.shape_values
+                              [k + (line * this->n_dofs_per_line()) +
+                               this->get_first_face_line_index(face_no)][i];
                       }
                   }
               }
diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h
index 4d6b8e6899..96fb938532 100644
--- a/include/deal.II/fe/fe_system.h
+++ b/include/deal.II/fe/fe_system.h
@@ -965,7 +965,7 @@ public:
    */
   virtual std::vector<std::pair<unsigned int, unsigned int>>
   hp_quad_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
-                         const unsigned int = 0) const override;
+                         const unsigned int face_no = 0) const override;
 
   /**
    * @copydoc FiniteElement::compare_for_domination()
@@ -1162,7 +1162,8 @@ private:
    */
   template <int structdim>
   std::vector<std::pair<unsigned int, unsigned int>>
-  hp_object_dof_identities(const FiniteElement<dim, spacedim> &fe_other) const;
+  hp_object_dof_identities(const FiniteElement<dim, spacedim> &fe_other,
+                           const unsigned int face_no = 0) const;
 
   /**
    * Usually: Fields of cell-independent data.
diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h
index 929d536ce4..a3c60daa3e 100644
--- a/include/deal.II/fe/fe_tools.h
+++ b/include/deal.II/fe/fe_tools.h
@@ -1246,7 +1246,8 @@ namespace FETools
       std::vector<std::pair<unsigned int, unsigned int>>
         &                                 face_system_to_component_table,
       const FiniteElement<dim, spacedim> &finite_element,
-      const bool                          do_tensor_product = true);
+      const bool                          do_tensor_product = true,
+      const unsigned int                  face_no           = 0 /*TODO*/);
 
   } // namespace Compositing
 
diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h
index 34c7d92789..cf8ac7d153 100644
--- a/include/deal.II/fe/fe_tools.templates.h
+++ b/include/deal.II/fe/fe_tools.templates.h
@@ -106,12 +106,17 @@ namespace FETools
       for (unsigned int i = 0; i < fes.size(); i++)
         if (multiplicities[i] > 0)
           {
+            // TODO: the implementation makes the assumption that all faces have
+            // the same number of dofs -> don't construct DPO but
+            // PrecomputedFiniteElementData
+            AssertDimension(fes[i]->n_unique_quads(), 1);
+
             multiplied_dofs_per_vertex +=
               fes[i]->n_dofs_per_vertex() * multiplicities[i];
             multiplied_dofs_per_line +=
               fes[i]->n_dofs_per_line() * multiplicities[i];
             multiplied_dofs_per_quad +=
-              fes[i]->n_dofs_per_quad() * multiplicities[i];
+              fes[i]->n_dofs_per_quad(0) * multiplicities[i];
             multiplied_dofs_per_hex +=
               fes[i]->n_dofs_per_hex() * multiplicities[i];
 
@@ -287,12 +292,11 @@ namespace FETools
           for (unsigned int base = 0; base < fes.size(); ++base)
             for (unsigned int m = 0; m < multiplicities[base]; ++m)
               for (unsigned int local_index = 0;
-                   local_index < fes[base]->n_dofs_per_quad();
+                   local_index < fes[base]->n_dofs_per_quad(quad_number);
                    ++local_index, ++total_index)
                 {
                   const unsigned int index_in_base =
-                    (fes[base]->n_dofs_per_quad() * quad_number + local_index +
-                     fes[base]->get_first_quad_index());
+                    local_index + fes[base]->get_first_quad_index(quad_number);
 
                   Assert(index_in_base < fes[base]->n_dofs_per_cell(),
                          ExcInternalError());
@@ -509,12 +513,11 @@ namespace FETools
                               comp_start +=
                               fes[base]->n_components() * do_tensor_product)
               for (unsigned int local_index = 0;
-                   local_index < fes[base]->n_dofs_per_quad();
+                   local_index < fes[base]->n_dofs_per_quad(quad_number);
                    ++local_index, ++total_index)
                 {
                   const unsigned int index_in_base =
-                    (fes[base]->n_dofs_per_quad() * quad_number + local_index +
-                     fes[base]->get_first_quad_index());
+                    local_index + fes[base]->get_first_quad_index(quad_number);
 
                   Assert(comp_start + fes[base]->n_components() <=
                            retval[total_index].size(),
@@ -773,13 +776,13 @@ namespace FETools
                               fe.base_element(base).n_components() *
                               do_tensor_product)
               for (unsigned int local_index = 0;
-                   local_index < fe.base_element(base).n_dofs_per_quad();
+                   local_index <
+                   fe.base_element(base).n_dofs_per_quad(quad_number);
                    ++local_index, ++total_index)
                 {
                   const unsigned int index_in_base =
-                    (fe.base_element(base).n_dofs_per_quad() * quad_number +
-                     local_index +
-                     fe.base_element(base).get_first_quad_index());
+                    local_index +
+                    fe.base_element(base).get_first_quad_index(quad_number);
 
                   system_to_base_table[total_index] =
                     std::make_pair(std::make_pair(base, m), index_in_base);
@@ -855,7 +858,8 @@ namespace FETools
       std::vector<std::pair<unsigned int, unsigned int>>
         &                                 face_system_to_component_table,
       const FiniteElement<dim, spacedim> &fe,
-      const bool                          do_tensor_product)
+      const bool                          do_tensor_product,
+      const unsigned int                  face_no)
     {
       // Initialize index tables. do this in the same way as done for the cell
       // tables, except that we now loop over the objects of faces
@@ -868,7 +872,8 @@ namespace FETools
       unsigned int total_index = 0;
       for (unsigned int vertex_number = 0;
            vertex_number <
-           ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0)
+           ReferenceCell::internal::Info::get_face(fe.reference_cell_type(),
+                                                   face_no)
              .n_vertices();
            ++vertex_number)
         {
@@ -906,12 +911,14 @@ namespace FETools
                     {
                       const unsigned int comp_in_base =
                         fe.base_element(base)
-                          .face_system_to_component_index(face_index_in_base)
+                          .face_system_to_component_index(face_index_in_base,
+                                                          face_no)
                           .first;
                       const unsigned int comp = comp_start + comp_in_base;
                       const unsigned int face_index_in_comp =
                         fe.base_element(base)
-                          .face_system_to_component_index(face_index_in_base)
+                          .face_system_to_component_index(face_index_in_base,
+                                                          face_no)
                           .second;
                       face_system_to_component_table[total_index] =
                         std::make_pair(comp, face_index_in_comp);
@@ -925,7 +932,8 @@ namespace FETools
       // 2. Lines
       for (unsigned int line_number = 0;
            line_number <
-           ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0)
+           ReferenceCell::internal::Info::get_face(fe.reference_cell_type(),
+                                                   face_no)
              .n_lines();
            ++line_number)
         {
@@ -947,7 +955,7 @@ namespace FETools
                      fe.base_element(base).get_first_line_index());
 
                   const unsigned int face_index_in_base =
-                    (fe.base_element(base).get_first_face_line_index() +
+                    (fe.base_element(base).get_first_face_line_index(face_no) +
                      fe.base_element(base).n_dofs_per_line() * line_number +
                      local_index);
 
@@ -958,12 +966,14 @@ namespace FETools
                     {
                       const unsigned int comp_in_base =
                         fe.base_element(base)
-                          .face_system_to_component_index(face_index_in_base)
+                          .face_system_to_component_index(face_index_in_base,
+                                                          face_no)
                           .first;
                       const unsigned int comp = comp_start + comp_in_base;
                       const unsigned int face_index_in_comp =
                         fe.base_element(base)
-                          .face_system_to_component_index(face_index_in_base)
+                          .face_system_to_component_index(face_index_in_base,
+                                                          face_no)
                           .second;
                       face_system_to_component_table[total_index] =
                         std::make_pair(comp, face_index_in_comp);
@@ -985,16 +995,16 @@ namespace FETools
                               fe.base_element(base).n_components() *
                               do_tensor_product)
               for (unsigned int local_index = 0;
-                   local_index < fe.base_element(base).n_dofs_per_quad();
+                   local_index < fe.base_element(base).n_dofs_per_quad(face_no);
                    ++local_index, ++total_index)
                 {
                   // do everything alike for this type of object
                   const unsigned int index_in_base =
                     (local_index +
-                     fe.base_element(base).get_first_quad_index());
+                     fe.base_element(base).get_first_quad_index(face_no));
 
                   const unsigned int face_index_in_base =
-                    (fe.base_element(base).get_first_face_quad_index() +
+                    (fe.base_element(base).get_first_face_quad_index(face_no) +
                      local_index);
 
                   face_system_to_base_table[total_index] =
@@ -1004,12 +1014,14 @@ namespace FETools
                     {
                       const unsigned int comp_in_base =
                         fe.base_element(base)
-                          .face_system_to_component_index(face_index_in_base)
+                          .face_system_to_component_index(face_index_in_base,
+                                                          face_no)
                           .first;
                       const unsigned int comp = comp_start + comp_in_base;
                       const unsigned int face_index_in_comp =
                         fe.base_element(base)
-                          .face_system_to_component_index(face_index_in_base)
+                          .face_system_to_component_index(face_index_in_base,
+                                                          face_no)
                           .second;
                       face_system_to_component_table[total_index] =
                         std::make_pair(comp, face_index_in_comp);
@@ -1019,7 +1031,7 @@ namespace FETools
                       non_primitive_index;
                 }
         }
-      Assert(total_index == fe.n_dofs_per_face(), ExcInternalError());
+      Assert(total_index == fe.n_dofs_per_face(face_no), ExcInternalError());
       Assert(total_index == face_system_to_component_table.size(),
              ExcInternalError());
       Assert(total_index == face_system_to_base_table.size(),
@@ -1934,11 +1946,13 @@ namespace FETools
     const unsigned int face_fine,
     const double       threshold)
   {
+    const unsigned int face_no = face_coarse;
+
     Assert(face_coarse == 0, ExcNotImplemented());
     Assert(face_fine == 0, ExcNotImplemented());
 
     const unsigned int nc     = GeometryInfo<dim>::max_children_per_face;
-    const unsigned int n      = fe.n_dofs_per_face();
+    const unsigned int n      = fe.n_dofs_per_face(face_no);
     const unsigned int nd     = fe.n_components();
     const unsigned int degree = fe.degree;
 
@@ -1961,9 +1975,9 @@ namespace FETools
     {
       unsigned int face_dof = 0;
       for (unsigned int i = 0;
-           i <
-           ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0)
-             .n_vertices();
+           i < ReferenceCell::internal::Info::get_face(fe.reference_cell_type(),
+                                                       face_no)
+                 .n_vertices();
            ++i)
         {
           const unsigned int offset_c =
@@ -1980,10 +1994,9 @@ namespace FETools
             }
         }
 
-      for (unsigned int i = 1;
-           i <=
-           ReferenceCell::internal::Info::get_face(fe.reference_cell_type(), 0)
-             .n_lines();
+      for (unsigned int i = 1; i <= ReferenceCell::internal::Info::get_face(
+                                      fe.reference_cell_type(), face_no)
+                                      .n_lines();
            ++i)
         {
           const unsigned int offset_c =
@@ -2004,18 +2017,16 @@ namespace FETools
 
       if (dim == 3)
         {
-          const unsigned int offset_c =
-            fe.get_first_quad_index() + face_coarse * fe.n_dofs_per_quad();
-          const unsigned int offset_f =
-            fe.get_first_quad_index() + face_fine * fe.n_dofs_per_quad();
-          for (unsigned int j = 0; j < fe.n_dofs_per_quad(); ++j)
+          const unsigned int offset_c = fe.get_first_quad_index(face_coarse);
+          const unsigned int offset_f = fe.get_first_quad_index(face_fine);
+          for (unsigned int j = 0; j < fe.n_dofs_per_quad(face_no); ++j)
             {
               face_c_dofs[face_dof] = offset_c + j;
               face_f_dofs[face_dof] = offset_f + j;
               ++face_dof;
             }
         }
-      Assert(face_dof == fe.n_dofs_per_face(), ExcInternalError());
+      Assert(face_dof == fe.n_dofs_per_face(face_no), ExcInternalError());
     }
 
     // Set up meshes, one with a single
diff --git a/include/deal.II/fe/fe_tools_extrapolate.templates.h b/include/deal.II/fe/fe_tools_extrapolate.templates.h
index 9f622655ad..b85649a2bf 100644
--- a/include/deal.II/fe/fe_tools_extrapolate.templates.h
+++ b/include/deal.II/fe/fe_tools_extrapolate.templates.h
@@ -1412,9 +1412,8 @@ namespace FETools
       compute_all_non_local_data(dof2, u2_relevant);
 
       // exclude dofs on more refined ghosted cells
-      const FiniteElement<dim, spacedim> &fe            = dof2.get_fe();
-      const unsigned int                  dofs_per_face = fe.n_dofs_per_face();
-      if (dofs_per_face > 0)
+      const FiniteElement<dim, spacedim> &fe = dof2.get_fe();
+      if (fe.max_dofs_per_face() > 0)
         {
           const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
           std::vector<types::global_dof_index> indices(dofs_per_cell);
@@ -1432,7 +1431,8 @@ namespace FETools
                       const typename DoFHandler<dim, spacedim>::cell_iterator
                         neighbor = cell->neighbor(face);
                       if (neighbor->level() != cell->level())
-                        for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        for (unsigned int i = 0; i < fe.n_dofs_per_face(face);
+                             ++i)
                           {
                             const types::global_dof_index index =
                               indices[fe.face_to_cell_index(i, face)];
diff --git a/source/fe/fe.cc b/source/fe/fe.cc
index 53a850375a..abd2bd4655 100644
--- a/source/fe/fe.cc
+++ b/source/fe/fe.cc
@@ -107,20 +107,27 @@ FiniteElement<dim, spacedim>::FiniteElement(
       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[0][j] =
-          std::pair<unsigned, unsigned>(0, j);
+      face_system_to_component_table.resize(this->n_unique_faces());
+      for (unsigned int f = 0; f < this->n_unique_faces(); ++f)
+        {
+          face_system_to_component_table[f].resize(this->n_dofs_per_face(f));
+          for (unsigned int j = 0; j < this->n_dofs_per_face(f); ++j)
+            face_system_to_component_table[f][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[0][j] = std::make_pair(std::make_pair(0U, 0U), j);
+  face_system_to_base_table.resize(this->n_unique_faces());
+  for (unsigned int f = 0; f < this->n_unique_faces(); ++f)
+    {
+      face_system_to_base_table[f].resize(this->n_dofs_per_face(f));
+      for (unsigned int j = 0; j < this->n_dofs_per_face(f); ++j)
+        face_system_to_base_table[f][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);
@@ -141,17 +148,26 @@ FiniteElement<dim, spacedim>::FiniteElement(
                                   FullMatrix<double>());
     }
 
+
   if (dim == 3)
     {
-      adjust_quad_dof_index_for_face_orientation_table.resize(1);
+      adjust_quad_dof_index_for_face_orientation_table.resize(
+        this->n_unique_quads());
 
-      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);
+      for (unsigned int f = 0; f < this->n_unique_quads(); ++f)
+        {
+          adjust_quad_dof_index_for_face_orientation_table[f] = Table<2, int>(
+            this->n_dofs_per_quad(f),
+            ReferenceCell::internal::Info::get_cell(this->reference_cell_type())
+                  .face_reference_cell_type(f) == ReferenceCell::Type::Quad ?
+              8 :
+              6);
+          adjust_quad_dof_index_for_face_orientation_table[f].fill(0);
+        }
     }
 
-  unit_face_support_points.resize(1);
-  generalized_face_support_points.resize(1);
+  unit_face_support_points.resize(this->n_unique_faces());
+  generalized_face_support_points.resize(this->n_unique_faces());
 }
 
 
@@ -558,7 +574,7 @@ FiniteElement<dim, spacedim>::face_to_cell_index(const unsigned int face_index,
   const auto &refence_cell =
     ReferenceCell::internal::Info::get_cell(this->reference_cell_type());
 
-  AssertIndexRange(face_index, this->n_dofs_per_face());
+  AssertIndexRange(face_index, this->n_dofs_per_face(face));
   AssertIndexRange(face, refence_cell.n_faces());
 
   // TODO: we could presumably solve the 3d case below using the
@@ -576,7 +592,7 @@ FiniteElement<dim, spacedim>::face_to_cell_index(const unsigned int face_index,
   // other than standard orientation
   if ((face_orientation != true) || (face_flip != false) ||
       (face_rotation != false))
-    Assert((this->n_dofs_per_line() <= 1) && (this->n_dofs_per_quad() <= 1),
+    Assert((this->n_dofs_per_line() <= 1) && (this->n_dofs_per_quad(face) <= 1),
            ExcMessage(
              "The function in this base class can not handle this case. "
              "Rather, the derived class you are using must provide "
@@ -585,7 +601,7 @@ FiniteElement<dim, spacedim>::face_to_cell_index(const unsigned int face_index,
 
   // we need to distinguish between DoFs on vertices, lines and in 3d quads.
   // do so in a sequence of if-else statements
-  if (face_index < this->get_first_face_line_index())
+  if (face_index < this->get_first_face_line_index(face))
     // DoF is on a vertex
     {
       // get the number of the vertex on the face that corresponds to this DoF,
@@ -604,12 +620,13 @@ FiniteElement<dim, spacedim>::face_to_cell_index(const unsigned int face_index,
                 this->n_dofs_per_vertex() +
               dof_index_on_vertex);
     }
-  else if (face_index < this->get_first_face_quad_index())
+  else if (face_index < this->get_first_face_quad_index(face))
     // DoF is on a face
     {
       // do the same kind of translation as before. we need to only consider
       // DoFs on the lines, i.e., ignoring those on the vertices
-      const unsigned int index = face_index - this->get_first_face_line_index();
+      const unsigned int index =
+        face_index - this->get_first_face_line_index(face);
 
       const unsigned int face_line         = index / this->n_dofs_per_line();
       const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
@@ -629,10 +646,10 @@ FiniteElement<dim, spacedim>::face_to_cell_index(const unsigned int face_index,
       Assert(dim >= 3, ExcInternalError());
 
       // ignore vertex and line dofs
-      const unsigned int index = face_index - this->get_first_face_quad_index();
+      const unsigned int index =
+        face_index - this->get_first_face_quad_index(face);
 
-      return (this->get_first_quad_index() + face * this->n_dofs_per_quad() +
-              index);
+      return (this->get_first_quad_index(face) + index);
     }
 }
 
@@ -642,13 +659,11 @@ template <int dim, int spacedim>
 unsigned int
 FiniteElement<dim, spacedim>::adjust_quad_dof_index_for_face_orientation(
   const unsigned int index,
-  const unsigned int face_no,
+  const unsigned int face,
   const bool         face_orientation,
   const bool         face_flip,
   const bool         face_rotation) const
 {
-  (void)face_no;
-
   // general template for 1D and 2D: not
   // implemented. in fact, the function
   // shouldn't even be called unless we are
@@ -666,12 +681,20 @@ FiniteElement<dim, spacedim>::adjust_quad_dof_index_for_face_orientation(
   // in 3d), so we don't need the table, but
   // the function should also not have been
   // called
-  AssertIndexRange(index, this->n_dofs_per_quad());
-  Assert(adjust_quad_dof_index_for_face_orientation_table[0].n_elements() ==
-           8 * this->n_dofs_per_quad(),
+  AssertIndexRange(index, this->n_dofs_per_quad(face));
+  Assert(adjust_quad_dof_index_for_face_orientation_table
+             [this->n_unique_quads() == 1 ? 0 : face]
+               .n_elements() ==
+           (ReferenceCell::internal::Info::get_cell(this->reference_cell_type())
+                  .face_reference_cell_type(face) == ReferenceCell::Type::Quad ?
+              8 :
+              6) *
+             this->n_dofs_per_quad(face),
          ExcInternalError());
-  return index + adjust_quad_dof_index_for_face_orientation_table[0](
-                   index, 4 * face_orientation + 2 * face_flip + face_rotation);
+  return index +
+         adjust_quad_dof_index_for_face_orientation_table
+           [this->n_unique_quads() == 1 ? 0 : face](
+             index, 4 * face_orientation + 2 * face_flip + face_rotation);
 }
 
 
@@ -818,8 +841,14 @@ bool
 FiniteElement<dim, spacedim>::constraints_are_implemented(
   const internal::SubfaceCase<dim> &subface_case) const
 {
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   if (subface_case == internal::SubfaceCase<dim>::case_isotropic)
-    return (this->n_dofs_per_face() == 0) || (interface_constraints.m() != 0);
+    return (this->n_dofs_per_face(face_no) == 0) ||
+           (interface_constraints.m() != 0);
   else
     return false;
 }
@@ -840,6 +869,12 @@ const FullMatrix<double> &
 FiniteElement<dim, spacedim>::constraints(
   const internal::SubfaceCase<dim> &subface_case) const
 {
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+  (void)face_no;
+
   (void)subface_case;
   Assert(subface_case == internal::SubfaceCase<dim>::case_isotropic,
          ExcMessage("Constraints for this element are only implemented "
@@ -847,7 +882,8 @@ FiniteElement<dim, spacedim>::constraints(
                     "(which is always the case in 2d, and in 3d requires "
                     "that the neighboring cell of a coarse cell presents "
                     "exactly four children on the common face)."));
-  Assert((this->n_dofs_per_face() == 0) || (interface_constraints.m() != 0),
+  Assert((this->n_dofs_per_face(face_no) == 0) ||
+           (interface_constraints.m() != 0),
          ExcMessage("The finite element for which you try to obtain "
                     "hanging node constraints does not appear to "
                     "implement them."));
@@ -866,17 +902,22 @@ template <int dim, int spacedim>
 TableIndices<2>
 FiniteElement<dim, spacedim>::interface_constraints_size() const
 {
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   switch (dim)
     {
       case 1:
         return {0U, 0U};
       case 2:
         return {this->n_dofs_per_vertex() + 2 * this->n_dofs_per_line(),
-                this->n_dofs_per_face()};
+                this->n_dofs_per_face(face_no)};
       case 3:
         return {5 * this->n_dofs_per_vertex() + 12 * this->n_dofs_per_line() +
-                  4 * this->n_dofs_per_quad(),
-                this->n_dofs_per_face()};
+                  4 * this->n_dofs_per_quad(face_no),
+                this->n_dofs_per_face(face_no)};
       default:
         Assert(false, ExcNotImplemented());
     }
@@ -1069,38 +1110,45 @@ FiniteElement<dim, spacedim>::unit_support_point(const unsigned int index) const
 template <int dim, int spacedim>
 const std::vector<Point<dim - 1>> &
 FiniteElement<dim, spacedim>::get_unit_face_support_points(
-  const unsigned int) const
+  const unsigned int face_no) const
 {
   // a finite element may define
   // support points, but only if
   // there are as many as there are
   // degrees of freedom on a face
-  Assert((unit_face_support_points[0].size() == 0) ||
-           (unit_face_support_points[0].size() == this->n_dofs_per_face()),
+  Assert((unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+            .size() == 0) ||
+           (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+              .size() == this->n_dofs_per_face(face_no)),
          ExcInternalError());
-  return unit_face_support_points[0];
+  return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no];
 }
 
 
 
 template <int dim, int spacedim>
 bool
-FiniteElement<dim, spacedim>::has_face_support_points(const unsigned int) const
+FiniteElement<dim, spacedim>::has_face_support_points(
+  const unsigned int face_no) const
 {
-  return (unit_face_support_points[0].size() != 0);
+  return (unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+            .size() != 0);
 }
 
 
 
 template <int dim, int spacedim>
 Point<dim - 1>
-FiniteElement<dim, spacedim>::unit_face_support_point(const unsigned int index,
-                                                      const unsigned int) const
+FiniteElement<dim, spacedim>::unit_face_support_point(
+  const unsigned int index,
+  const unsigned int face_no) const
 {
-  AssertIndexRange(index, this->n_dofs_per_face());
-  Assert(unit_face_support_points[0].size() == this->n_dofs_per_face(),
+  AssertIndexRange(index, this->n_dofs_per_face(face_no));
+  Assert(unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+             .size() == this->n_dofs_per_face(face_no),
          ExcFEHasNoSupportPoints());
-  return unit_face_support_points[0][index];
+  return unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+                                 [index];
 }
 
 
diff --git a/source/fe/fe_abf.cc b/source/fe/fe_abf.cc
index c55552234b..0b11c719cf 100644
--- a/source/fe/fe_abf.cc
+++ b/source/fe/fe_abf.cc
@@ -88,17 +88,24 @@ FE_ABF<dim>::FE_ABF(const unsigned int deg)
 
   initialize_restriction();
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   // TODO[TL]: for anisotropic refinement we will probably need a table of
   // submatrices with an array for each refine case
   std::vector<FullMatrix<double>> face_embeddings(
     1 << (dim - 1),
-    FullMatrix<double>(this->n_dofs_per_face(), this->n_dofs_per_face()));
+    FullMatrix<double>(this->n_dofs_per_face(face_no),
+                       this->n_dofs_per_face(face_no)));
   // TODO: Something goes wrong there. The error of the least squares fit
   // is to large ...
   // FETools::compute_face_embedding_matrices(*this, face_embeddings.data(), 0,
   // 0);
-  this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(),
-                                     this->n_dofs_per_face());
+  this->interface_constraints.reinit((1 << (dim - 1)) *
+                                       this->n_dofs_per_face(face_no),
+                                     this->n_dofs_per_face(face_no));
   unsigned int target_row = 0;
   for (const auto &face_embedding : face_embeddings)
     for (unsigned int i = 0; i < face_embedding.m(); ++i)
@@ -153,6 +160,11 @@ FE_ABF<dim>::initialize_support_points(const unsigned int deg)
   QGauss<dim>        cell_quadrature(deg + 2);
   const unsigned int n_interior_points = cell_quadrature.size();
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   unsigned int n_face_points = (dim > 1) ? 1 : 0;
   // compute (deg+1)^(dim-1)
   for (unsigned int d = 1; d < dim; ++d)
@@ -160,7 +172,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[0].resize(n_face_points);
+  this->generalized_face_support_points[face_no].resize(n_face_points);
 
 
   // These might be required when the faces contribution is computed
@@ -194,7 +206,8 @@ 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[0][k] = face_points.point(k);
+          this->generalized_face_support_points[face_no][k] =
+            face_points.point(k);
           // Compute its quadrature
           // contribution for each
           // moment.
@@ -381,7 +394,8 @@ FE_ABF<dim>::initialize_restriction()
           for (unsigned int k = 0; k < n_face_points; ++k)
             for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
                  ++i_child)
-              for (unsigned int i_face = 0; i_face < this->n_dofs_per_face();
+              for (unsigned int i_face = 0;
+                   i_face < this->n_dofs_per_face(face);
                    ++i_face)
                 {
                   // The quadrature
@@ -389,13 +403,12 @@ FE_ABF<dim>::initialize_restriction()
                   // subcell are NOT
                   // transformed, so we
                   // have to do it here.
-                  this->restriction[iso][child](face * this->n_dofs_per_face() +
-                                                  i_face,
-                                                i_child) +=
+                  this->restriction[iso][child](
+                    face * this->n_dofs_per_face(face) + i_face, i_child) +=
                     Utilities::fixed_power<dim - 1>(.5) * q_sub.weight(k) *
                     cached_values_face(i_child, k) *
                     this->shape_value_component(
-                      face * this->n_dofs_per_face() + i_face,
+                      face * this->n_dofs_per_face(face) + i_face,
                       q_sub.point(k),
                       GeometryInfo<dim>::unit_normal_direction[face]);
                 }
@@ -419,9 +432,14 @@ FE_ABF<dim>::initialize_restriction()
       polynomials[dd] = new AnisotropicPolynomials<dim>(poly);
     }
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   QGauss<dim>        q_cell(rt_order + 1);
   const unsigned int start_cell_dofs =
-    GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face();
+    GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
 
   // Store shape values, since the
   // evaluation suffers if not
@@ -561,14 +579,19 @@ FE_ABF<dim>::convert_generalized_support_point_values_to_dof_values(
     for (unsigned int k = 0; k < n_face_points; ++k)
       for (unsigned int i = 0; i < boundary_weights.size(1); ++i)
         {
-          nodal_values[i + face * this->n_dofs_per_face()] +=
+          nodal_values[i + face * this->n_dofs_per_face(face)] +=
             boundary_weights(k, i) *
             support_point_values[face * n_face_points + k][GeometryInfo<
               dim>::unit_normal_direction[face]];
         }
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   const unsigned int start_cell_dofs =
-    GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face();
+    GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
   const unsigned int start_cell_points =
     GeometryInfo<dim>::faces_per_cell * n_face_points;
 
diff --git a/source/fe/fe_bdm.cc b/source/fe/fe_bdm.cc
index dd1e4b469e..2efea8a775 100644
--- a/source/fe/fe_bdm.cc
+++ b/source/fe/fe_bdm.cc
@@ -80,12 +80,17 @@ FE_BDM<dim>::FE_BDM(const unsigned int deg)
   this->reinit_restriction_and_prolongation_matrices(true, true);
   FETools::compute_embedding_matrices(*this, this->prolongation, true, 1.);
 
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
   for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
-    face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face());
+    face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
+                              this->n_dofs_per_face(face_no));
   FETools::compute_face_embedding_matrices(*this, face_embeddings, 0, 0, 1.);
-  this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(),
-                                     this->n_dofs_per_face());
+  this->interface_constraints.reinit((1 << (dim - 1)) *
+                                       this->n_dofs_per_face(face_no),
+                                     this->n_dofs_per_face(face_no));
   unsigned int target_row = 0;
   for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
     for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
@@ -155,15 +160,15 @@ FE_BDM<dim>::convert_generalized_support_point_values_to_dof_values(
       // initialize_support_points()
       if (test_values_face.size() == 0)
         {
-          for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+          for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
             nodal_values[dbase + i] =
               support_point_values[pbase + i]
                                   [GeometryInfo<dim>::unit_normal_direction[f]];
-          pbase += this->n_dofs_per_face();
+          pbase += this->n_dofs_per_face(f);
         }
       else
         {
-          for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+          for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
             {
               double s = 0.;
               for (unsigned int k = 0; k < test_values_face.size(); ++k)
@@ -175,11 +180,16 @@ FE_BDM<dim>::convert_generalized_support_point_values_to_dof_values(
             }
           pbase += test_values_face.size();
         }
-      dbase += this->n_dofs_per_face();
+      dbase += this->n_dofs_per_face(f);
     }
 
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+  (void)face_no;
+
   AssertDimension(dbase,
-                  this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell);
+                  this->n_dofs_per_face(face_no) *
+                    GeometryInfo<dim>::faces_per_cell);
   AssertDimension(pbase,
                   this->generalized_support_points.size() -
                     test_values_cell.size());
@@ -332,10 +342,15 @@ FE_BDM<dim>::initialize_support_points(const unsigned int deg)
   // considered later. In 2D, we can use point values.
   QGauss<dim - 1> face_points(deg + 1);
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   // Copy the quadrature formula to the face points.
-  this->generalized_face_support_points[0].resize(face_points.size());
+  this->generalized_face_support_points[face_no].resize(face_points.size());
   for (unsigned int k = 0; k < face_points.size(); ++k)
-    this->generalized_face_support_points[0][k] = face_points.point(k);
+    this->generalized_face_support_points[face_no][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
@@ -356,12 +371,13 @@ FE_BDM<dim>::initialize_support_points(const unsigned int deg)
        k < face_points.size() * GeometryInfo<dim>::faces_per_cell;
        ++k)
     this->generalized_support_points[k] = faces.point(
-      k + QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
-                                                   0,
-                                                   true,
-                                                   false,
-                                                   false,
-                                                   this->n_dofs_per_face()));
+      k +
+      QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
+                                               0,
+                                               true,
+                                               false,
+                                               false,
+                                               this->n_dofs_per_face(face_no)));
 
   // Currently, for backward compatibility, we do not use moments, but
   // point values on faces in 2D. In 3D, this is impossible, since the
diff --git a/source/fe/fe_bernstein.cc b/source/fe/fe_bernstein.cc
index 8dee487809..86525d8659 100644
--- a/source/fe/fe_bernstein.cc
+++ b/source/fe/fe_bernstein.cc
@@ -112,11 +112,11 @@ FE_Bernstein<dim, spacedim>::get_subface_interpolation_matrix(
   const FiniteElement<dim, spacedim> &x_source_fe,
   const unsigned int                  subface,
   FullMatrix<double> &                interpolation_matrix,
-  const unsigned int) const
+  const unsigned int                  face_no) const
 {
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
 
   // see if source is a Bernstein element
   if (const FE_Bernstein<dim, spacedim> *source_fe =
@@ -124,9 +124,9 @@ FE_Bernstein<dim, spacedim>::get_subface_interpolation_matrix(
     {
       // have this test in here since a table of size 2x0 reports its size as
       // 0x0
-      Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+      Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
              ExcDimensionMismatch(interpolation_matrix.n(),
-                                  this->n_dofs_per_face()));
+                                  this->n_dofs_per_face(face_no)));
 
       // Make sure that the element for which the DoFs should be constrained
       // is the one with the higher polynomial degree.  Actually the procedure
@@ -134,13 +134,13 @@ FE_Bernstein<dim, spacedim>::get_subface_interpolation_matrix(
       // produced in that case might lead to problems in the hp procedures,
       // which use this method.
       Assert(
-        this->n_dofs_per_face() <= source_fe->n_dofs_per_face(),
+        this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
         (typename FiniteElement<dim,
                                 spacedim>::ExcInterpolationNotImplemented()));
 
       const Quadrature<dim - 1> quad_face_support(
         FE_Q<dim, spacedim>(QIterated<1>(QTrapez<1>(), source_fe->degree))
-          .get_unit_face_support_points());
+          .get_unit_face_support_points(face_no));
 
       // Rule of thumb for FP accuracy, that can be expected for a given
       // polynomial degree.  This value is used to cut off values close to
@@ -163,10 +163,10 @@ FE_Bernstein<dim, spacedim>::get_subface_interpolation_matrix(
                                               0,
                                               subface);
 
-      for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i)
+      for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
         {
           const Point<dim> &p = subface_quadrature.point(i);
-          for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
+          for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
             {
               double matrix_entry =
                 this->shape_value(this->face_to_cell_index(j, 0), p);
@@ -185,11 +185,11 @@ FE_Bernstein<dim, spacedim>::get_subface_interpolation_matrix(
 
       // make sure that the row sum of each of the matrices is 1 at this
       // point. this must be so since the shape functions sum up to 1
-      for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j)
+      for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
         {
           double sum = 0.;
 
-          for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+          for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
             sum += interpolation_matrix(j, i);
 
           Assert(std::fabs(sum - 1) < eps, ExcInternalError());
@@ -236,7 +236,7 @@ FE_Bernstein<dim, spacedim>::hp_vertex_dof_identities(
       // equivalencies to be recorded
       return std::vector<std::pair<unsigned int, unsigned int>>();
     }
-  else if (fe_other.n_dofs_per_face() == 0)
+  else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
     {
       // if the other element has no elements on faces at all,
       // then it would be impossible to enforce any kind of
diff --git a/source/fe/fe_enriched.cc b/source/fe/fe_enriched.cc
index 6f1cc3d2c5..f470b2e008 100644
--- a/source/fe/fe_enriched.cc
+++ b/source/fe/fe_enriched.cc
@@ -440,7 +440,6 @@ 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[0].resize(this->n_dofs_per_face());
 
     FETools::Compositing::build_cell_tables(this->system_to_base_table,
                                             this->system_to_component_table,
@@ -448,11 +447,21 @@ FE_Enriched<dim, spacedim>::initialize(
                                             *this,
                                             false);
 
-    FETools::Compositing::build_face_tables(
-      this->face_system_to_base_table[0],
-      this->face_system_to_component_table[0],
-      *this,
-      false);
+    this->face_system_to_component_table.resize(this->n_unique_faces());
+
+    for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no)
+      {
+        this->face_system_to_component_table[0].resize(
+          this->n_dofs_per_face(face_no));
+
+
+        FETools::Compositing::build_face_tables(
+          this->face_system_to_base_table[face_no],
+          this->face_system_to_component_table[face_no],
+          *this,
+          false,
+          face_no);
+      }
   }
 
   // restriction and prolongation matrices are built on demand
@@ -468,8 +477,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[0] = fe_system->unit_face_support_points[0];
+    this->unit_support_points      = fe_system->unit_support_points;
+    this->unit_face_support_points = fe_system->unit_face_support_points;
   }
 
   // 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 b9534a584e..d673916566 100644
--- a/source/fe/fe_face.cc
+++ b/source/fe/fe_face.cc
@@ -163,16 +163,16 @@ FE_FaceQ<dim, spacedim>::get_subface_interpolation_matrix(
   const FiniteElement<dim, spacedim> &x_source_fe,
   const unsigned int                  subface,
   FullMatrix<double> &                interpolation_matrix,
-  const unsigned int) const
+  const unsigned int                  face_no) const
 {
   // this function is similar to the respective method in FE_Q
 
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+                              this->n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
 
   // see if source is a FaceQ element
   if (const FE_FaceQ<dim, spacedim> *source_fe =
@@ -184,13 +184,13 @@ FE_FaceQ<dim, spacedim>::get_subface_interpolation_matrix(
       // produced in that case might lead to problems in the hp procedures,
       // which use this method.
       Assert(
-        this->n_dofs_per_face() <= source_fe->n_dofs_per_face(),
+        this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
         (typename FiniteElement<dim,
                                 spacedim>::ExcInterpolationNotImplemented()));
 
       // generate a quadrature with the unit face support points.
       const Quadrature<dim - 1> face_quadrature(
-        source_fe->get_unit_face_support_points());
+        source_fe->get_unit_face_support_points(face_no));
 
       // Rule of thumb for FP accuracy, that can be expected for a given
       // polynomial degree.  This value is used to cut off values close to
@@ -199,7 +199,7 @@ FE_FaceQ<dim, spacedim>::get_subface_interpolation_matrix(
 
       // compute the interpolation matrix by simply taking the value at the
       // support points.
-      for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i)
+      for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
         {
           const Point<dim - 1> p =
             subface == numbers::invalid_unsigned_int ?
@@ -207,7 +207,7 @@ FE_FaceQ<dim, spacedim>::get_subface_interpolation_matrix(
               GeometryInfo<dim - 1>::child_to_cell_coordinates(
                 face_quadrature.point(i), subface);
 
-          for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
+          for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
             {
               double matrix_entry = this->poly_space.compute_value(j, p);
 
@@ -225,11 +225,11 @@ FE_FaceQ<dim, spacedim>::get_subface_interpolation_matrix(
 
       // make sure that the row sum of each of the matrices is 1 at this
       // point. this must be so since the shape functions sum up to 1
-      for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j)
+      for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
         {
           double sum = 0.;
 
-          for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+          for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
             sum += interpolation_matrix(j, i);
 
           Assert(std::fabs(sum - 1) < eps, ExcInternalError());
@@ -254,7 +254,7 @@ FE_FaceQ<dim, spacedim>::has_support_on_face(
   const unsigned int shape_index,
   const unsigned int face_index) const
 {
-  return (face_index == (shape_index / this->n_dofs_per_face()));
+  return (face_index == (shape_index / this->n_dofs_per_face(face_index)));
 }
 
 
@@ -340,7 +340,8 @@ FE_FaceQ<dim, spacedim>::hp_line_dof_identities(
           // equivalencies to be recorded
           return std::vector<std::pair<unsigned int, unsigned int>>();
         }
-      else if (fe_other.n_dofs_per_face() == 0)
+      else if (fe_other.n_unique_faces() == 1 &&
+               fe_other.n_dofs_per_face(0) == 0)
         {
           // if the other element has no elements on faces at all,
           // then it would be impossible to enforce any kind of
@@ -418,7 +419,8 @@ FE_FaceQ<dim, spacedim>::hp_quad_dof_identities(
           // equivalencies to be recorded
           return std::vector<std::pair<unsigned int, unsigned int>>();
         }
-      else if (fe_other.n_dofs_per_face() == 0)
+      else if (fe_other.n_unique_faces() == 1 &&
+               fe_other.n_dofs_per_face(0) == 0)
         {
           // if the other element has no elements on faces at all,
           // then it would be impossible to enforce any kind of
@@ -574,15 +576,17 @@ FE_FaceQ<1, spacedim>::get_subface_interpolation_matrix(
   const FiniteElement<1, spacedim> &x_source_fe,
   const unsigned int /*subface*/,
   FullMatrix<double> &interpolation_matrix,
-  const unsigned int) const
+  const unsigned int  face_no) const
 {
   (void)x_source_fe;
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+  (void)face_no;
+
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+                              this->n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
   interpolation_matrix(0, 0) = 1.;
 }
 
@@ -798,7 +802,7 @@ FE_FaceP<dim, spacedim>::has_support_on_face(
   const unsigned int shape_index,
   const unsigned int face_index) const
 {
-  return (face_index == (shape_index / this->n_dofs_per_face()));
+  return (face_index == (shape_index / this->n_dofs_per_face(face_index)));
 }
 
 
@@ -888,16 +892,16 @@ FE_FaceP<dim, spacedim>::get_subface_interpolation_matrix(
   const FiniteElement<dim, spacedim> &x_source_fe,
   const unsigned int                  subface,
   FullMatrix<double> &                interpolation_matrix,
-  const unsigned int) const
+  const unsigned int                  face_no) const
 {
   // this function is similar to the respective method in FE_Q
 
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+                              this->n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
 
   // see if source is a FaceP element
   if (const FE_FaceP<dim, spacedim> *source_fe =
@@ -909,7 +913,7 @@ FE_FaceP<dim, spacedim>::get_subface_interpolation_matrix(
       // produced in that case might lead to problems in the hp procedures,
       // which use this method.
       Assert(
-        this->n_dofs_per_face() <= source_fe->n_dofs_per_face(),
+        this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
         (typename FiniteElement<dim,
                                 spacedim>::ExcInterpolationNotImplemented()));
 
@@ -924,7 +928,7 @@ FE_FaceP<dim, spacedim>::get_subface_interpolation_matrix(
       const double eps = 2e-13 * (this->degree + 1) * (dim - 1);
 
       FullMatrix<double> mass(face_quadrature.size(),
-                              source_fe->n_dofs_per_face());
+                              source_fe->n_dofs_per_face(face_no));
 
       for (unsigned int k = 0; k < face_quadrature.size(); ++k)
         {
@@ -934,18 +938,18 @@ FE_FaceP<dim, spacedim>::get_subface_interpolation_matrix(
               GeometryInfo<dim - 1>::child_to_cell_coordinates(
                 face_quadrature.point(k), subface);
 
-          for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j)
+          for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
             mass(k, j) = source_fe->poly_space.compute_value(j, p);
         }
 
       Householder<double> H(mass);
       Vector<double>      v_in(face_quadrature.size());
-      Vector<double>      v_out(source_fe->n_dofs_per_face());
+      Vector<double>      v_out(source_fe->n_dofs_per_face(face_no));
 
 
       // compute the interpolation matrix by evaluating on the fine side and
       // then solving the least squares problem
-      for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+      for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
         {
           for (unsigned int k = 0; k < face_quadrature.size(); ++k)
             {
@@ -960,7 +964,7 @@ FE_FaceP<dim, spacedim>::get_subface_interpolation_matrix(
           (void)result;
           Assert(result < 1e-12, FETools::ExcLeastSquaresError(result));
 
-          for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j)
+          for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
             {
               double matrix_entry = v_out(j);
 
@@ -995,7 +999,7 @@ FE_FaceP<dim, spacedim>::get_constant_modes() const
 {
   Table<2, bool> constant_modes(1, this->n_dofs_per_cell());
   for (unsigned int face : GeometryInfo<dim>::face_indices())
-    constant_modes(0, face * this->n_dofs_per_face()) = true;
+    constant_modes(0, face * this->n_dofs_per_face(face)) = true;
   return std::pair<Table<2, bool>, std::vector<unsigned int>>(
     constant_modes, std::vector<unsigned int>(1, 0));
 }
diff --git a/source/fe/fe_nedelec.cc b/source/fe/fe_nedelec.cc
index 6f46551ccb..caa4f1f39d 100644
--- a/source/fe/fe_nedelec.cc
+++ b/source/fe/fe_nedelec.cc
@@ -108,8 +108,14 @@ FE_Nedelec<dim>::FE_Nedelec(const unsigned int order)
 #endif
   FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
-    face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face());
+    face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
+                              this->n_dofs_per_face(face_no));
 
   FETools::compute_face_embedding_matrices<dim, double>(
     *this,
@@ -128,14 +134,15 @@ FE_Nedelec<dim>::FE_Nedelec(const unsigned int order)
 
       case 2:
         {
-          this->interface_constraints.reinit(2 * this->n_dofs_per_face(),
-                                             this->n_dofs_per_face());
+          this->interface_constraints.reinit(2 * this->n_dofs_per_face(face_no),
+                                             this->n_dofs_per_face(face_no));
 
           for (unsigned int i = 0; i < GeometryInfo<2>::max_children_per_face;
                ++i)
-            for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
-              for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
-                this->interface_constraints(i * this->n_dofs_per_face() + j,
+            for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
+              for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
+                this->interface_constraints(i * this->n_dofs_per_face(face_no) +
+                                              j,
                                             k) = face_embeddings[i](j, k);
 
           break;
@@ -143,16 +150,16 @@ FE_Nedelec<dim>::FE_Nedelec(const unsigned int order)
 
       case 3:
         {
-          this->interface_constraints.reinit(4 * (this->n_dofs_per_face() -
-                                                  this->degree),
-                                             this->n_dofs_per_face());
+          this->interface_constraints.reinit(
+            4 * (this->n_dofs_per_face(face_no) - this->degree),
+            this->n_dofs_per_face(face_no));
 
           unsigned int target_row = 0;
 
           for (unsigned int i = 0; i < 2; ++i)
             for (unsigned int j = this->degree; j < 2 * this->degree;
                  ++j, ++target_row)
-              for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
+              for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
                 this->interface_constraints(target_row, k) =
                   face_embeddings[2 * i](j, k);
 
@@ -160,7 +167,7 @@ FE_Nedelec<dim>::FE_Nedelec(const unsigned int order)
             for (unsigned int j = 3 * this->degree;
                  j < GeometryInfo<3>::lines_per_face * this->degree;
                  ++j, ++target_row)
-              for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
+              for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
                 this->interface_constraints(target_row, k) =
                   face_embeddings[i](j, k);
 
@@ -169,7 +176,8 @@ FE_Nedelec<dim>::FE_Nedelec(const unsigned int order)
               for (unsigned int k = i * this->degree;
                    k < (i + 1) * this->degree;
                    ++k, ++target_row)
-                for (unsigned int l = 0; l < this->n_dofs_per_face(); ++l)
+                for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
+                     ++l)
                   this->interface_constraints(target_row, l) =
                     face_embeddings[i + 2 * j](k, l);
 
@@ -178,7 +186,8 @@ FE_Nedelec<dim>::FE_Nedelec(const unsigned int order)
               for (unsigned int k = (i + 2) * this->degree;
                    k < (i + 3) * this->degree;
                    ++k, ++target_row)
-                for (unsigned int l = 0; l < this->n_dofs_per_face(); ++l)
+                for (unsigned int l = 0; l < this->n_dofs_per_face(face_no);
+                     ++l)
                   this->interface_constraints(target_row, l) =
                     face_embeddings[2 * i + j](k, l);
 
@@ -186,9 +195,9 @@ FE_Nedelec<dim>::FE_Nedelec(const unsigned int order)
                ++i)
             for (unsigned int j =
                    GeometryInfo<3>::lines_per_face * this->degree;
-                 j < this->n_dofs_per_face();
+                 j < this->n_dofs_per_face(face_no);
                  ++j, ++target_row)
-              for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
+              for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
                 this->interface_constraints(target_row, k) =
                   face_embeddings[i](j, k);
 
@@ -254,6 +263,11 @@ FE_Nedelec<2>::initialize_support_points(const unsigned int order)
 {
   const int dim = 2;
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   // Create polynomial basis.
   const std::vector<Polynomials::Polynomial<double>> &lobatto_polynomials =
     Polynomials::Lobatto::generate_complete_basis(order + 1);
@@ -273,11 +287,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[0].resize(n_edge_points);
+  this->generalized_face_support_points[face_no].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[0][q_point] =
+    this->generalized_face_support_points[face_no][q_point] =
       reference_edge_quadrature.point(q_point);
 
   if (order > 0)
@@ -310,7 +324,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[0][q_point](0));
+                this->generalized_face_support_points[face_no][q_point](0));
         }
 
       for (unsigned int q_point = 0; q_point < n_interior_points; ++q_point)
@@ -347,6 +361,11 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order)
 {
   const int dim = 3;
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   // Create polynomial basis.
   const std::vector<Polynomials::Polynomial<double>> &lobatto_polynomials =
     Polynomials::Lobatto::generate_complete_basis(order + 1);
@@ -380,8 +399,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[0].resize(4 * n_edge_points +
-                                                      n_face_points);
+      this->generalized_face_support_points[face_no].resize(4 * n_edge_points +
+                                                            n_face_points);
       this->generalized_support_points.resize(n_boundary_points +
                                               n_interior_points);
 
@@ -391,8 +410,9 @@ 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[0][line * n_edge_points +
-                                                     q_point] =
+            this
+              ->generalized_face_support_points[face_no][line * n_edge_points +
+                                                         q_point] =
               edge_quadrature.point(
                 QProjector<dim - 1>::DataSetDescriptor::face(
                   ReferenceCell::get_hypercube(dim - 1),
@@ -421,14 +441,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[0][q_point](1));
+                this->generalized_face_support_points[face_no][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[0][q_point + 4 * n_edge_points] =
+          this->generalized_face_support_points[face_no]
+                                               [q_point + 4 * n_edge_points] =
             reference_face_quadrature.point(q_point);
 
           for (unsigned int i = 0; i <= order; ++i)
@@ -437,24 +457,20 @@ 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[0][q_point +
-                                                             4 * n_edge_points](
-                      0)) *
+                    this->generalized_face_support_points
+                      [face_no][q_point + 4 * n_edge_points](0)) *
                   lobatto_polynomials[j + 2].value(
-                    this->generalized_face_support_points[0][q_point +
-                                                             4 * n_edge_points](
-                      1));
+                    this->generalized_face_support_points
+                      [face_no][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[0][q_point +
-                                                             4 * n_edge_points](
-                      1)) *
+                    this->generalized_face_support_points
+                      [face_no][q_point + 4 * n_edge_points](1)) *
                   lobatto_polynomials[j + 2].value(
-                    this->generalized_face_support_points[0][q_point +
-                                                             4 * n_edge_points](
-                      0));
+                    this->generalized_face_support_points
+                      [face_no][q_point + 4 * n_edge_points](0));
               }
         }
 
@@ -486,7 +502,7 @@ FE_Nedelec<3>::initialize_support_points(const unsigned int order)
 
   else
     {
-      this->generalized_face_support_points[0].resize(4 * n_edge_points);
+      this->generalized_face_support_points[face_no].resize(4 * n_edge_points);
       this->generalized_support_points.resize(
         GeometryInfo<dim>::lines_per_cell * n_edge_points);
 
@@ -495,8 +511,9 @@ 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[0][line * n_edge_points +
-                                                     q_point] =
+            this
+              ->generalized_face_support_points[face_no][line * n_edge_points +
+                                                         q_point] =
               edge_quadrature.point(
                 QProjector<dim - 1>::DataSetDescriptor::face(
                   ReferenceCell::get_hypercube(dim - 1),
@@ -2453,20 +2470,21 @@ void
 FE_Nedelec<dim>::get_face_interpolation_matrix(
   const FiniteElement<dim> &source,
   FullMatrix<double> &      interpolation_matrix,
-  const unsigned int) const
+  const unsigned int        face_no) const
 {
+  (void)face_no;
   // this is only implemented, if the
   // source FE is also a
   // Nedelec element
   AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
                 (dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
               (typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
-  Assert(interpolation_matrix.m() == source.n_dofs_per_face(),
+  Assert(interpolation_matrix.m() == source.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              source.n_dofs_per_face()));
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+                              source.n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
+                              this->n_dofs_per_face(face_no)));
 
   // ok, source is a Nedelec element, so
   // we will be able to do the work
@@ -2484,7 +2502,7 @@ FE_Nedelec<dim>::get_face_interpolation_matrix(
   // lead to problems in the
   // hp procedures, which use this
   // method.
-  Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
          (typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
   interpolation_matrix = 0;
 
@@ -2560,7 +2578,7 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
   const FiniteElement<dim> &source,
   const unsigned int        subface,
   FullMatrix<double> &      interpolation_matrix,
-  const unsigned int) const
+  const unsigned int        face_no) const
 {
   // this is only implemented, if the
   // source FE is also a
@@ -2568,12 +2586,12 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
   AssertThrow((source.get_name().find("FE_Nedelec<") == 0) ||
                 (dynamic_cast<const FE_Nedelec<dim> *>(&source) != nullptr),
               typename FiniteElement<dim>::ExcInterpolationNotImplemented());
-  Assert(interpolation_matrix.m() == source.n_dofs_per_face(),
+  Assert(interpolation_matrix.m() == source.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              source.n_dofs_per_face()));
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+                              source.n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
+                              this->n_dofs_per_face(face_no)));
 
   // ok, source is a Nedelec element, so
   // we will be able to do the work
@@ -2591,7 +2609,7 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
   // lead to problems in the
   // hp procedures, which use this
   // method.
-  Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
          (typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
   interpolation_matrix = 0.0;
   // Perform projection-based interpolation
@@ -2605,7 +2623,8 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
     {
       case 2:
         {
-          for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof)
+          for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
+               ++dof)
             for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
                  ++q_point)
               {
@@ -2653,7 +2672,8 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
               Vector<double> solution(source_fe.degree - 1);
               Vector<double> system_rhs(source_fe.degree - 1);
 
-              for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof)
+              for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
+                   ++dof)
                 {
                   system_rhs = 0.0;
 
@@ -2700,7 +2720,8 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
                                        {0.0, 1.0},
                                        {1.0, 1.0}};
 
-          for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof)
+          for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
+               ++dof)
             for (unsigned int q_point = 0; q_point < n_edge_quadrature_points;
                  ++q_point)
               {
@@ -2769,7 +2790,8 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
                                             GeometryInfo<dim>::lines_per_face);
               Vector<double>     tmp(GeometryInfo<dim>::lines_per_face);
 
-              for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof)
+              for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
+                   ++dof)
                 {
                   system_rhs = 0.0;
 
@@ -2892,7 +2914,8 @@ FE_Nedelec<dim>::get_subface_interpolation_matrix(
               system_rhs.reinit(system_matrix_inv.m(), 2);
               tmp.reinit(2);
 
-              for (unsigned int dof = 0; dof < this->n_dofs_per_face(); ++dof)
+              for (unsigned int dof = 0; dof < this->n_dofs_per_face(face_no);
+                   ++dof)
                 {
                   system_rhs = 0.0;
 
@@ -3112,6 +3135,11 @@ FE_Nedelec<dim>::convert_generalized_support_point_values_to_dof_values(
   const std::vector<Vector<double>> &support_point_values,
   std::vector<double> &              nodal_values) const
 {
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   const unsigned int deg = this->degree - 1;
   Assert(support_point_values.size() == this->generalized_support_points.size(),
          ExcDimensionMismatch(support_point_values.size(),
@@ -3186,7 +3214,8 @@ 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[0][q_point](0));
+                        this->generalized_face_support_points[face_no][q_point](
+                          0));
 
               FullMatrix<double> system_matrix_inv(this->degree - 1,
                                                    this->degree - 1);
@@ -3465,7 +3494,8 @@ 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[0][q_point](1));
+                        this->generalized_face_support_points[face_no][q_point](
+                          1));
 
               FullMatrix<double> system_matrix_inv(this->degree - 1,
                                                    this->degree - 1);
@@ -3544,10 +3574,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
-                              [0][q_point + 4 * n_edge_points](0)) *
+                              [face_no][q_point + 4 * n_edge_points](0)) *
                           lobatto_polynomials[j + 2].value(
                             this->generalized_face_support_points
-                              [0][q_point + 4 * n_edge_points](1));
+                              [face_no][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_nedelec_sz.cc b/source/fe/fe_nedelec_sz.cc
index 7a9d9bb646..f2a6e5d049 100644
--- a/source/fe/fe_nedelec_sz.cc
+++ b/source/fe/fe_nedelec_sz.cc
@@ -147,7 +147,9 @@ FE_NedelecSZ<dim, spacedim>::get_data(
   const unsigned int faces_per_cell    = GeometryInfo<dim>::faces_per_cell;
 
   const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell;
-  const unsigned int n_face_dofs = this->n_dofs_per_quad() * faces_per_cell;
+
+  // we assume that all quads have the same numer of dofs
+  const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell;
 
   const UpdateFlags  flags(data.update_each);
   const unsigned int n_q_points = quadrature.size();
@@ -1695,7 +1697,8 @@ FE_NedelecSZ<dim, spacedim>::fill_face_values(
           // Loop through quad points:
           for (unsigned int m = 0; m < faces_per_cell; ++m)
             {
-              const unsigned int shift_m(m * this->n_dofs_per_quad());
+              // we assume that all quads have the same numer of dofs
+              const unsigned int shift_m(m * this->n_dofs_per_quad(0));
               // Calculate the offsets for each face-based shape function:
               //
               // Type-1 (gradients)
diff --git a/source/fe/fe_poly_tensor.cc b/source/fe/fe_poly_tensor.cc
index d6f1794bf3..6df444121b 100644
--- a/source/fe/fe_poly_tensor.cc
+++ b/source/fe/fe_poly_tensor.cc
@@ -83,11 +83,11 @@ namespace internal
                 const unsigned int nn = cell->neighbor_face_no(f);
 
                 if (nn < GeometryInfo<dim>::faces_per_cell / 2)
-                  for (unsigned int j = 0; j < fe.n_dofs_per_face(); ++j)
+                  for (unsigned int j = 0; j < fe.n_dofs_per_face(f); ++j)
                     {
                       const unsigned int cell_j = fe.face_to_cell_index(j, f);
 
-                      Assert(f * fe.n_dofs_per_face() + j < face_sign.size(),
+                      Assert(f * fe.n_dofs_per_face(f) + j < face_sign.size(),
                              ExcInternalError());
                       Assert(mapping_kind.size() == 1 ||
                                cell_j < mapping_kind.size(),
@@ -98,7 +98,7 @@ namespace internal
                       if ((mapping_kind.size() > 1 ?
                              mapping_kind[cell_j] :
                              mapping_kind[0]) == mapping_raviart_thomas)
-                        face_sign[f * fe.n_dofs_per_face() + j] = -1.0;
+                        face_sign[f * fe.n_dofs_per_face(f) + j] = -1.0;
                     }
               }
           }
diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc
index 38c301c2ba..3d6cf65dd4 100644
--- a/source/fe/fe_q_base.cc
+++ b/source/fe/fe_q_base.cc
@@ -624,11 +624,11 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_subface_interpolation_matrix(
   const FiniteElement<dim, spacedim> &x_source_fe,
   const unsigned int                  subface,
   FullMatrix<double> &                interpolation_matrix,
-  const unsigned int) const
+  const unsigned int                  face_no) const
 {
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
 
   // see if source is a Q element
   if (const FE_Q_Base<PolynomialType, dim, spacedim> *source_fe =
@@ -637,9 +637,9 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_subface_interpolation_matrix(
     {
       // have this test in here since a table of size 2x0 reports its size as
       // 0x0
-      Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+      Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
              ExcDimensionMismatch(interpolation_matrix.n(),
-                                  this->n_dofs_per_face()));
+                                  this->n_dofs_per_face(face_no)));
 
       // Make sure that the element for which the DoFs should be constrained
       // is the one with the higher polynomial degree.  Actually the procedure
@@ -647,13 +647,13 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_subface_interpolation_matrix(
       // produced in that case might lead to problems in the hp procedures,
       // which use this method.
       Assert(
-        this->n_dofs_per_face() <= source_fe->n_dofs_per_face(),
+        this->n_dofs_per_face(face_no) <= source_fe->n_dofs_per_face(face_no),
         (typename FiniteElement<dim,
                                 spacedim>::ExcInterpolationNotImplemented()));
 
       // generate a point on this cell and evaluate the shape functions there
       const Quadrature<dim - 1> quad_face_support(
-        source_fe->get_unit_face_support_points());
+        source_fe->get_unit_face_support_points(face_no));
 
       // Rule of thumb for FP accuracy, that can be expected for a given
       // polynomial degree.  This value is used to cut off values close to
@@ -674,11 +674,11 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_subface_interpolation_matrix(
                                               quad_face_support,
                                               0,
                                               subface);
-      for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i)
+      for (unsigned int i = 0; i < source_fe->n_dofs_per_face(face_no); ++i)
         {
           const Point<dim> &p = subface_quadrature.point(i);
 
-          for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
+          for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
             {
               double matrix_entry =
                 this->shape_value(this->face_to_cell_index(j, 0), p);
@@ -697,11 +697,11 @@ FE_Q_Base<PolynomialType, dim, spacedim>::get_subface_interpolation_matrix(
 
       // make sure that the row sum of each of the matrices is 1 at this
       // point. this must be so since the shape functions sum up to 1
-      for (unsigned int j = 0; j < source_fe->n_dofs_per_face(); ++j)
+      for (unsigned int j = 0; j < source_fe->n_dofs_per_face(face_no); ++j)
         {
           double sum = 0.;
 
-          for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+          for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
             sum += interpolation_matrix(j, i);
 
           Assert(std::fabs(sum - 1) < eps, ExcInternalError());
@@ -750,7 +750,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::hp_vertex_dof_identities(
       // equivalencies to be recorded
       return std::vector<std::pair<unsigned int, unsigned int>>();
     }
-  else if (fe_other.n_dofs_per_face() == 0)
+  else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
     {
       // if the other element has no elements on faces at all,
       // then it would be impossible to enforce any kind of
@@ -815,7 +815,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::hp_line_dof_identities(
       // equivalencies to be recorded
       return std::vector<std::pair<unsigned int, unsigned int>>();
     }
-  else if (fe_other.n_dofs_per_face() == 0)
+  else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
     {
       // if the other element has no elements on faces at all,
       // then it would be impossible to enforce any kind of
@@ -886,7 +886,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::hp_quad_dof_identities(
       // equivalencies to be recorded
       return std::vector<std::pair<unsigned int, unsigned int>>();
     }
-  else if (fe_other.n_dofs_per_face() == 0)
+  else if (fe_other.n_unique_faces() == 1 && fe_other.n_dofs_per_face(0) == 0)
     {
       // if the other element has no elements on faces at all,
       // then it would be impossible to enforce any kind of
@@ -946,7 +946,12 @@ FE_Q_Base<PolynomialType, dim, spacedim>::initialize_unit_face_support_points(
   if (dim == 1)
     return;
 
-  this->unit_face_support_points[0].resize(
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
+  this->unit_face_support_points[face_no].resize(
     Utilities::fixed_power<dim - 1>(q_degree + 1));
 
   // find renumbering of faces and assign from values of quadrature
@@ -962,9 +967,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[0].resize(support_quadrature.size());
+  this->unit_face_support_points[face_no].resize(support_quadrature.size());
   for (unsigned int k = 0; k < support_quadrature.size(); ++k)
-    this->unit_face_support_points[0][face_index_map[k]] =
+    this->unit_face_support_points[face_no][face_index_map[k]] =
       support_quadrature.point(k);
 }
 
@@ -979,12 +984,17 @@ FE_Q_Base<PolynomialType, dim, spacedim>::
   if (dim < 3)
     return;
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   Assert(this->adjust_quad_dof_index_for_face_orientation_table[0]
-             .n_elements() == 8 * this->n_dofs_per_quad(),
+             .n_elements() == 8 * this->n_dofs_per_quad(face_no),
          ExcInternalError());
 
   const unsigned int n = q_degree - 1;
-  Assert(n * n == this->n_dofs_per_quad(), ExcInternalError());
+  Assert(n * n == this->n_dofs_per_quad(face_no), ExcInternalError());
 
   // the dofs on a face are connected to a n x n matrix. for example, for
   // degree==4 we have the following dofs on a quad
@@ -1004,34 +1014,42 @@ FE_Q_Base<PolynomialType, dim, spacedim>::
   // rotated and mirrored numbers.
 
 
-  for (unsigned int local = 0; local < this->n_dofs_per_quad(); ++local)
+  for (unsigned int local = 0; local < this->n_dofs_per_quad(face_no); ++local)
     // face support points are in lexicographic ordering with x running
     // fastest. invert that (y running fastest)
     {
       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[0](local, 0) =
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+                                                                      0) =
         j + i * n - local;
       // face_orientation=false, face_flip=false, face_rotation=true
-      this->adjust_quad_dof_index_for_face_orientation_table[0](local, 1) =
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](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[0](local, 2) =
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](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[0](local, 3) =
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](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[0](local, 4) = 0;
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+                                                                      4) = 0;
       // face_orientation=true,  face_flip=false, face_rotation=true
-      this->adjust_quad_dof_index_for_face_orientation_table[0](local, 5) =
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](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[0](local, 6) =
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](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[0](local, 7) =
+      this->adjust_quad_dof_index_for_face_orientation_table[face_no](local,
+                                                                      7) =
         (n - 1 - j) + i * n - local;
     }
 
@@ -1052,7 +1070,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::face_to_cell_index(
   const bool         face_flip,
   const bool         face_rotation) const
 {
-  AssertIndexRange(face_index, this->n_dofs_per_face());
+  AssertIndexRange(face_index, this->n_dofs_per_face(face));
   AssertIndexRange(face, GeometryInfo<dim>::faces_per_cell);
 
   // TODO: we could presumably solve the 3d case below using the
@@ -1066,7 +1084,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::face_to_cell_index(
 
   // we need to distinguish between DoFs on vertices, lines and in 3d quads.
   // do so in a sequence of if-else statements
-  if (face_index < this->get_first_face_line_index())
+  if (face_index < this->get_first_face_line_index(face))
     // DoF is on a vertex
     {
       // get the number of the vertex on the face that corresponds to this DoF,
@@ -1082,12 +1100,13 @@ FE_Q_Base<PolynomialType, dim, spacedim>::face_to_cell_index(
                 this->n_dofs_per_vertex() +
               dof_index_on_vertex);
     }
-  else if (face_index < this->get_first_face_quad_index())
+  else if (face_index < this->get_first_face_quad_index(face))
     // DoF is on a face
     {
       // do the same kind of translation as before. we need to only consider
       // DoFs on the lines, i.e., ignoring those on the vertices
-      const unsigned int index = face_index - this->get_first_face_line_index();
+      const unsigned int index =
+        face_index - this->get_first_face_line_index(face);
 
       const unsigned int face_line         = index / this->n_dofs_per_line();
       const unsigned int dof_index_on_line = index % this->n_dofs_per_line();
@@ -1142,17 +1161,17 @@ FE_Q_Base<PolynomialType, dim, spacedim>::face_to_cell_index(
       Assert(dim >= 3, ExcInternalError());
 
       // ignore vertex and line dofs
-      const unsigned int index = face_index - this->get_first_face_quad_index();
+      const unsigned int index =
+        face_index - this->get_first_face_quad_index(face);
 
       // the same is true here as above for the 3d case -- someone will
       // just have to draw a bunch of pictures. in the meantime,
       // we can implement the Q2 case in which it is simple
-      Assert((this->n_dofs_per_quad() <= 1) ||
+      Assert((this->n_dofs_per_quad(face) <= 1) ||
                ((face_orientation == true) && (face_flip == false) &&
                 (face_rotation == false)),
              ExcNotImplemented());
-      return (this->get_first_quad_index() + face * this->n_dofs_per_quad() +
-              index);
+      return (this->get_first_quad_index(face) + index);
     }
 }
 
@@ -1544,7 +1563,8 @@ FE_Q_Base<PolynomialType, dim, spacedim>::has_support_on_face(
 
   // first, special-case interior shape functions, since they have no support
   // no-where on the boundary
-  if (((dim == 2) && (shape_index >= this->get_first_quad_index())) ||
+  if (((dim == 2) &&
+       (shape_index >= this->get_first_quad_index(0 /*first quad*/))) ||
       ((dim == 3) && (shape_index >= this->get_first_hex_index())))
     return false;
 
@@ -1565,7 +1585,7 @@ FE_Q_Base<PolynomialType, dim, spacedim>::has_support_on_face(
 
       return false;
     }
-  else if (shape_index < this->get_first_quad_index())
+  else if (shape_index < this->get_first_quad_index(0 /*first quad*/))
     // ok, dof is on a line
     {
       const unsigned int line_index =
@@ -1596,7 +1616,8 @@ FE_Q_Base<PolynomialType, dim, spacedim>::has_support_on_face(
     // dof is on a quad
     {
       const unsigned int quad_index =
-        (shape_index - this->get_first_quad_index()) / this->n_dofs_per_quad();
+        (shape_index - this->get_first_quad_index(0)) /
+        this->n_dofs_per_quad(face_index); // this won't work
       Assert(static_cast<signed int>(quad_index) <
                static_cast<signed int>(GeometryInfo<dim>::quads_per_cell),
              ExcInternalError());
diff --git a/source/fe/fe_q_hierarchical.cc b/source/fe/fe_q_hierarchical.cc
index 3af55b0c28..670ec26d5e 100644
--- a/source/fe/fe_q_hierarchical.cc
+++ b/source/fe/fe_q_hierarchical.cc
@@ -307,8 +307,13 @@ FE_Q_Hierarchical<dim>::hp_quad_dof_identities(
   // one is an FE_Nothing.
   if (dynamic_cast<const FE_Q_Hierarchical<dim> *>(&fe_other) != nullptr)
     {
-      const unsigned int this_dpq  = this->n_dofs_per_quad();
-      const unsigned int other_dpq = fe_other.n_dofs_per_quad();
+      // TODO: the implementation makes the assumption that all faces have the
+      // same number of dofs
+      AssertDimension(this->n_unique_faces(), 1);
+      const unsigned int face_no = 0;
+
+      const unsigned int this_dpq  = this->n_dofs_per_quad(face_no);
+      const unsigned int other_dpq = fe_other.n_dofs_per_quad(face_no);
 
       // we deal with hierarchical 1d polynomials where dofs are enumerated
       // increasingly. Thus we return a vector of pairs for the first N-1, where
@@ -927,7 +932,7 @@ void
 FE_Q_Hierarchical<dim>::get_face_interpolation_matrix(
   const FiniteElement<dim> &x_source_fe,
   FullMatrix<double> &      interpolation_matrix,
-  const unsigned int) const
+  const unsigned int        face_no) const
 {
   // this is only implemented, if the
   // source FE is also a
@@ -938,12 +943,12 @@ FE_Q_Hierarchical<dim>::get_face_interpolation_matrix(
                  nullptr),
               (typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
 
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+                              this->n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
 
   // ok, source is a Q_Hierarchical element, so
   // we will be able to do the work
@@ -962,7 +967,7 @@ FE_Q_Hierarchical<dim>::get_face_interpolation_matrix(
   // lead to problems in the
   // hp procedures, which use this
   // method.
-  Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
          (typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
   interpolation_matrix = 0;
 
@@ -976,7 +981,7 @@ FE_Q_Hierarchical<dim>::get_face_interpolation_matrix(
           // element, which corresponds to 1 on diagonal of the matrix.
           // DoFs which correspond to higher polynomials
           // are zeroed (zero rows in the matrix).
-          for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+          for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
             interpolation_matrix(i, i) = 1;
 
           break;
@@ -1016,7 +1021,7 @@ FE_Q_Hierarchical<dim>::get_subface_interpolation_matrix(
   const FiniteElement<dim> &x_source_fe,
   const unsigned int        subface,
   FullMatrix<double> &      interpolation_matrix,
-  const unsigned int) const
+  const unsigned int        face_no) const
 {
   // this is only implemented, if the
   // source FE is also a
@@ -1027,12 +1032,12 @@ FE_Q_Hierarchical<dim>::get_subface_interpolation_matrix(
                  nullptr),
               (typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
 
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+                              this->n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
 
   // ok, source is a Q_Hierarchical element, so
   // we will be able to do the work
@@ -1050,7 +1055,7 @@ FE_Q_Hierarchical<dim>::get_subface_interpolation_matrix(
   // lead to problems in the
   // hp procedures, which use this
   // method.
-  Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
          (typename FiniteElement<dim>::ExcInterpolationNotImplemented()));
 
   switch (dim)
@@ -1065,7 +1070,8 @@ FE_Q_Hierarchical<dim>::get_subface_interpolation_matrix(
                   interpolation_matrix(1, 0) = 0.5;
                   interpolation_matrix(1, 1) = 0.5;
 
-                  for (unsigned int dof = 2; dof < this->n_dofs_per_face();)
+                  for (unsigned int dof = 2;
+                       dof < this->n_dofs_per_face(face_no);)
                     {
                       interpolation_matrix(1, dof) = -1.0;
                       dof                          = dof + 2;
@@ -1073,14 +1079,16 @@ FE_Q_Hierarchical<dim>::get_subface_interpolation_matrix(
 
                   int factorial_i = 1;
 
-                  for (unsigned int i = 2; i < this->n_dofs_per_face(); ++i)
+                  for (unsigned int i = 2; i < this->n_dofs_per_face(face_no);
+                       ++i)
                     {
                       interpolation_matrix(i, i) = std::pow(0.5, i);
                       factorial_i *= i;
                       int factorial_j  = factorial_i;
                       int factorial_ij = 1;
 
-                      for (unsigned int j = i + 1; j < this->n_dofs_per_face();
+                      for (unsigned int j = i + 1;
+                           j < this->n_dofs_per_face(face_no);
                            ++j)
                         {
                           factorial_ij *= j - i;
@@ -1106,7 +1114,8 @@ FE_Q_Hierarchical<dim>::get_subface_interpolation_matrix(
                   interpolation_matrix(0, 0) = 0.5;
                   interpolation_matrix(0, 1) = 0.5;
 
-                  for (unsigned int dof = 2; dof < this->n_dofs_per_face();)
+                  for (unsigned int dof = 2;
+                       dof < this->n_dofs_per_face(face_no);)
                     {
                       interpolation_matrix(0, dof) = -1.0;
                       dof                          = dof + 2;
@@ -1116,14 +1125,16 @@ FE_Q_Hierarchical<dim>::get_subface_interpolation_matrix(
 
                   int factorial_i = 1;
 
-                  for (unsigned int i = 2; i < this->n_dofs_per_face(); ++i)
+                  for (unsigned int i = 2; i < this->n_dofs_per_face(face_no);
+                       ++i)
                     {
                       interpolation_matrix(i, i) = std::pow(0.5, i);
                       factorial_i *= i;
                       int factorial_j  = factorial_i;
                       int factorial_ij = 1;
 
-                      for (unsigned int j = i + 1; j < this->n_dofs_per_face();
+                      for (unsigned int j = i + 1;
+                           j < this->n_dofs_per_face(face_no);
                            ++j)
                         {
                           factorial_ij *= j - i;
@@ -1903,12 +1914,17 @@ FE_Q_Hierarchical<dim>::initialize_generalized_face_support_points()
 {
   const unsigned int codim = dim - 1;
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   // number of points: (degree+1)^codim
   unsigned int n = this->degree + 1;
   for (unsigned int i = 1; i < codim; ++i)
     n *= this->degree + 1;
 
-  this->generalized_face_support_points[0].resize(n);
+  this->generalized_face_support_points[face_no].resize(n);
 
   Point<codim> p;
 
@@ -1941,7 +1957,8 @@ FE_Q_Hierarchical<dim>::initialize_generalized_face_support_points()
               else
                 p(2) = .5;
             }
-          this->generalized_face_support_points[0][face_renumber[k++]] = p;
+          this->generalized_face_support_points[face_no][face_renumber[k++]] =
+            p;
         }
 }
 
@@ -2029,7 +2046,7 @@ FE_Q_Hierarchical<dim>::hierarchic_to_fe_q_hierarchical_numbering(
           for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
             h2l[next_index++] = n + 2 + i;
           // inside quad
-          Assert(fe.n_dofs_per_quad() ==
+          Assert(fe.n_dofs_per_quad(0 /*only one quad in 2D*/) ==
                    fe.n_dofs_per_line() * fe.n_dofs_per_line(),
                  ExcInternalError());
           for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
@@ -2086,8 +2103,14 @@ FE_Q_Hierarchical<dim>::hierarchic_to_fe_q_hierarchical_numbering(
           for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
             h2l[next_index++] = (2 + i) * n2 + n + 1;
 
+          // TODO: the implementation makes the assumption that all faces have
+          // the same number of dofs
+          AssertDimension(fe.n_unique_faces(), 1);
+          const unsigned int face_no = 0;
+          (void)face_no;
+
           // inside quads
-          Assert(fe.n_dofs_per_quad() ==
+          Assert(fe.n_dofs_per_quad(face_no) ==
                    fe.n_dofs_per_line() * fe.n_dofs_per_line(),
                  ExcInternalError());
           // left face
@@ -2117,7 +2140,7 @@ FE_Q_Hierarchical<dim>::hierarchic_to_fe_q_hierarchical_numbering(
 
           // inside hex
           Assert(fe.n_dofs_per_hex() ==
-                   fe.n_dofs_per_quad() * fe.n_dofs_per_line(),
+                   fe.n_dofs_per_quad(face_no) * fe.n_dofs_per_line(),
                  ExcInternalError());
           for (unsigned int i = 0; i < fe.n_dofs_per_line(); ++i)
             for (unsigned int j = 0; j < fe.n_dofs_per_line(); ++j)
@@ -2192,7 +2215,8 @@ FE_Q_Hierarchical<dim>::has_support_on_face(const unsigned int shape_index,
   // shape functions, since they
   // have no support no-where on
   // the boundary
-  if (((dim == 2) && (shape_index >= this->get_first_quad_index())) ||
+  if (((dim == 2) && (shape_index >=
+                      this->get_first_quad_index(0 /*only one quad in 2D*/))) ||
       ((dim == 3) && (shape_index >= this->get_first_hex_index())))
     return false;
 
@@ -2214,7 +2238,7 @@ FE_Q_Hierarchical<dim>::has_support_on_face(const unsigned int shape_index,
           return true;
       return false;
     }
-  else if (shape_index < this->get_first_quad_index())
+  else if (shape_index < this->get_first_quad_index(0))
     // ok, dof is on a line
     {
       const unsigned int line_index =
@@ -2231,7 +2255,8 @@ FE_Q_Hierarchical<dim>::has_support_on_face(const unsigned int shape_index,
     // dof is on a quad
     {
       const unsigned int quad_index =
-        (shape_index - this->get_first_quad_index()) / this->n_dofs_per_quad();
+        (shape_index - this->get_first_quad_index(0 /*first quad*/)) /
+        this->n_dofs_per_quad(face_index);
       Assert(static_cast<signed int>(quad_index) <
                static_cast<signed int>(GeometryInfo<dim>::quads_per_cell),
              ExcInternalError());
diff --git a/source/fe/fe_raviart_thomas.cc b/source/fe/fe_raviart_thomas.cc
index 5b4334ae6c..78d3271a9a 100644
--- a/source/fe/fe_raviart_thomas.cc
+++ b/source/fe/fe_raviart_thomas.cc
@@ -87,17 +87,24 @@ FE_RaviartThomas<dim>::FE_RaviartThomas(const unsigned int deg)
   FETools::compute_embedding_matrices(*this, this->prolongation);
   initialize_restriction();
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   // TODO[TL]: for anisotropic refinement we will probably need a table of
   // submatrices with an array for each refine case
   FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
   for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
-    face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face());
+    face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
+                              this->n_dofs_per_face(face_no));
   FETools::compute_face_embedding_matrices<dim, double>(*this,
                                                         face_embeddings,
                                                         0,
                                                         0);
-  this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(),
-                                     this->n_dofs_per_face());
+  this->interface_constraints.reinit((1 << (dim - 1)) *
+                                       this->n_dofs_per_face(face_no),
+                                     this->n_dofs_per_face(face_no));
   unsigned int target_row = 0;
   for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
     for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
@@ -152,6 +159,11 @@ FE_RaviartThomas<dim>::initialize_support_points(const unsigned int deg)
   QGauss<dim>        cell_quadrature(deg + 1);
   const unsigned int n_interior_points = (deg > 0) ? cell_quadrature.size() : 0;
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   unsigned int n_face_points = (dim > 1) ? 1 : 0;
   // compute (deg+1)^(dim-1)
   for (unsigned int d = 1; d < dim; ++d)
@@ -160,7 +172,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[0].resize(n_face_points);
+  this->generalized_face_support_points[face_no].resize(n_face_points);
 
   // Number of the point being entered
   unsigned int current = 0;
@@ -173,12 +185,10 @@ FE_RaviartThomas<dim>::initialize_support_points(const unsigned int deg)
 
       boundary_weights.reinit(n_face_points, legendre.n());
 
-      //       Assert (face_points.size() == this->n_dofs_per_face(),
-      //            ExcInternalError());
-
       for (unsigned int k = 0; k < n_face_points; ++k)
         {
-          this->generalized_face_support_points[0][k] = face_points.point(k);
+          this->generalized_face_support_points[face_no][k] =
+            face_points.point(k);
           // Compute its quadrature
           // contribution for each
           // moment.
@@ -323,7 +333,8 @@ FE_RaviartThomas<dim>::initialize_restriction()
           for (unsigned int k = 0; k < n_face_points; ++k)
             for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
                  ++i_child)
-              for (unsigned int i_face = 0; i_face < this->n_dofs_per_face();
+              for (unsigned int i_face = 0;
+                   i_face < this->n_dofs_per_face(face);
                    ++i_face)
                 {
                   // The quadrature
@@ -331,13 +342,12 @@ FE_RaviartThomas<dim>::initialize_restriction()
                   // subcell are NOT
                   // transformed, so we
                   // have to do it here.
-                  this->restriction[iso][child](face * this->n_dofs_per_face() +
-                                                  i_face,
-                                                i_child) +=
+                  this->restriction[iso][child](
+                    face * this->n_dofs_per_face(face) + i_face, i_child) +=
                     Utilities::fixed_power<dim - 1>(.5) * q_sub.weight(k) *
                     cached_values_on_face(i_child, k) *
                     this->shape_value_component(
-                      face * this->n_dofs_per_face() + i_face,
+                      face * this->n_dofs_per_face(face) + i_face,
                       q_sub.point(k),
                       GeometryInfo<dim>::unit_normal_direction[face]);
                 }
@@ -362,9 +372,14 @@ FE_RaviartThomas<dim>::initialize_restriction()
       polynomials[dd] = std::make_unique<AnisotropicPolynomials<dim>>(poly);
     }
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   QGauss<dim>        q_cell(this->degree);
   const unsigned int start_cell_dofs =
-    GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face();
+    GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
 
   // Store shape values, since the
   // evaluation suffers if not
@@ -512,14 +527,19 @@ FE_RaviartThomas<dim>::convert_generalized_support_point_values_to_dof_values(
     for (unsigned int k = 0; k < n_face_points; ++k)
       for (unsigned int i = 0; i < boundary_weights.size(1); ++i)
         {
-          nodal_values[i + face * this->n_dofs_per_face()] +=
+          nodal_values[i + face * this->n_dofs_per_face(face)] +=
             boundary_weights(k, i) *
             support_point_values[face * n_face_points + k](
               GeometryInfo<dim>::unit_normal_direction[face]);
         }
 
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   const unsigned int start_cell_dofs =
-    GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face();
+    GeometryInfo<dim>::faces_per_cell * this->n_dofs_per_face(face_no);
   const unsigned int start_cell_points =
     GeometryInfo<dim>::faces_per_cell * n_face_points;
 
diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc
index 90ca1195c6..4ca716558b 100644
--- a/source/fe/fe_raviart_thomas_nodal.cc
+++ b/source/fe/fe_raviart_thomas_nodal.cc
@@ -82,19 +82,27 @@ FE_RaviartThomasNodal<dim>::FE_RaviartThomasNodal(const unsigned int deg)
       for (unsigned int i = 0; i < nc; ++i)
         this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
     }
+
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   // Fill prolongation matrices with embedding operators
   FETools::compute_embedding_matrices(*this, this->prolongation);
   // TODO[TL]: for anisotropic refinement we will probably need a table of
   // submatrices with an array for each refine case
   FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
   for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
-    face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face());
+    face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
+                              this->n_dofs_per_face(face_no));
   FETools::compute_face_embedding_matrices<dim, double>(*this,
                                                         face_embeddings,
                                                         0,
                                                         0);
-  this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(),
-                                     this->n_dofs_per_face());
+  this->interface_constraints.reinit((1 << (dim - 1)) *
+                                       this->n_dofs_per_face(face_no),
+                                     this->n_dofs_per_face(face_no));
   unsigned int target_row = 0;
   for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
     for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
@@ -147,8 +155,14 @@ template <int dim>
 void
 FE_RaviartThomasNodal<dim>::initialize_support_points(const unsigned int deg)
 {
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   this->generalized_support_points.resize(this->n_dofs_per_cell());
-  this->generalized_face_support_points[0].resize(this->n_dofs_per_face());
+  this->generalized_face_support_points[face_no].resize(
+    this->n_dofs_per_face(face_no));
 
   // Number of the point being entered
   unsigned int current = 0;
@@ -162,25 +176,28 @@ FE_RaviartThomasNodal<dim>::initialize_support_points(const unsigned int deg)
   if (dim > 1)
     {
       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[0][k] = face_points.point(k);
+      Assert(face_points.size() == this->n_dofs_per_face(face_no),
+             ExcInternalError());
+      for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
+        this->generalized_face_support_points[face_no][k] =
+          face_points.point(k);
       Quadrature<dim> faces =
         QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
                                               face_points);
-      for (unsigned int k = 0;
-           k < this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
+      for (unsigned int k = 0; k < this->n_dofs_per_face(face_no) *
+                                     GeometryInfo<dim>::faces_per_cell;
            ++k)
-        this->generalized_support_points[k] = faces.point(
-          k +
-          QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
-                                                   0,
-                                                   true,
-                                                   false,
-                                                   false,
-                                                   this->n_dofs_per_face()));
-
-      current = this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
+        this->generalized_support_points[k] =
+          faces.point(k + QProjector<dim>::DataSetDescriptor::face(
+                            this->reference_cell_type(),
+                            0,
+                            true,
+                            false,
+                            false,
+                            this->n_dofs_per_face(face_no)));
+
+      current =
+        this->n_dofs_per_face(face_no) * GeometryInfo<dim>::faces_per_cell;
     }
 
   if (deg == 0)
@@ -329,9 +346,9 @@ FE_RaviartThomasNodal<dim>::
   unsigned int fbase = 0;
   unsigned int f     = 0;
   for (; f < GeometryInfo<dim>::faces_per_cell;
-       ++f, fbase += this->n_dofs_per_face())
+       ++f, fbase += this->n_dofs_per_face(f))
     {
-      for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+      for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
         {
           nodal_values[fbase + i] = support_point_values[fbase + i](
             GeometryInfo<dim>::unit_normal_direction[f]);
@@ -463,7 +480,7 @@ template <int dim>
 std::vector<std::pair<unsigned int, unsigned int>>
 FE_RaviartThomasNodal<dim>::hp_quad_dof_identities(
   const FiniteElement<dim> &fe_other,
-  const unsigned int) const
+  const unsigned int        face_no) const
 {
   // we can presently only compute
   // these identities if both FEs are
@@ -479,8 +496,10 @@ FE_RaviartThomasNodal<dim>::hp_quad_dof_identities(
 
       // this works exactly like the line
       // case above
-      const unsigned int p = this->n_dofs_per_quad();
-      const unsigned int q = fe_q_other->n_dofs_per_quad();
+      const unsigned int p = this->n_dofs_per_quad(face_no);
+
+      AssertDimension(fe_q_other->n_unique_faces(), 1);
+      const unsigned int q = fe_q_other->n_dofs_per_quad(0);
 
       std::vector<std::pair<unsigned int, unsigned int>> identities;
 
@@ -575,7 +594,7 @@ void
 FE_RaviartThomasNodal<dim>::get_face_interpolation_matrix(
   const FiniteElement<dim> &x_source_fe,
   FullMatrix<double> &      interpolation_matrix,
-  const unsigned int) const
+  const unsigned int        face_no) const
 {
   // this is only implemented, if the
   // source FE is also a
@@ -585,12 +604,12 @@ FE_RaviartThomasNodal<dim>::get_face_interpolation_matrix(
                    &x_source_fe) != nullptr),
               typename FiniteElement<dim>::ExcInterpolationNotImplemented());
 
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+                              this->n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
 
   // ok, source is a RaviartThomasNodal element, so
   // we will be able to do the work
@@ -608,7 +627,7 @@ FE_RaviartThomasNodal<dim>::get_face_interpolation_matrix(
   // lead to problems in the
   // hp procedures, which use this
   // method.
-  Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
          typename FiniteElement<dim>::ExcInterpolationNotImplemented());
 
   // generate a quadrature
@@ -618,7 +637,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[0]);
+    source_fe.generalized_face_support_points[face_no]);
 
   // Rule of thumb for FP accuracy,
   // that can be expected for a
@@ -635,11 +654,11 @@ FE_RaviartThomasNodal<dim>::get_face_interpolation_matrix(
                                      quad_face_support,
                                      0);
 
-  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(); ++i)
+  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
     {
       const Point<dim> &p = face_projection.point(i);
 
-      for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
+      for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
         {
           double matrix_entry =
             this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
@@ -664,11 +683,11 @@ FE_RaviartThomasNodal<dim>::get_face_interpolation_matrix(
   // this point. this must be so
   // since the shape functions sum up
   // to 1
-  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(); ++j)
+  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
     {
       double sum = 0.;
 
-      for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+      for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
         sum += interpolation_matrix(j, i);
 
       Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
@@ -683,7 +702,7 @@ FE_RaviartThomasNodal<dim>::get_subface_interpolation_matrix(
   const FiniteElement<dim> &x_source_fe,
   const unsigned int        subface,
   FullMatrix<double> &      interpolation_matrix,
-  const unsigned int) const
+  const unsigned int        face_no) const
 {
   // this is only implemented, if the
   // source FE is also a
@@ -693,12 +712,12 @@ FE_RaviartThomasNodal<dim>::get_subface_interpolation_matrix(
                    &x_source_fe) != nullptr),
               typename FiniteElement<dim>::ExcInterpolationNotImplemented());
 
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+                              this->n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
 
   // ok, source is a RaviartThomasNodal element, so
   // we will be able to do the work
@@ -716,7 +735,7 @@ FE_RaviartThomasNodal<dim>::get_subface_interpolation_matrix(
   // lead to problems in the
   // hp procedures, which use this
   // method.
-  Assert(this->n_dofs_per_face() <= source_fe.n_dofs_per_face(),
+  Assert(this->n_dofs_per_face(face_no) <= source_fe.n_dofs_per_face(face_no),
          typename FiniteElement<dim>::ExcInterpolationNotImplemented());
 
   // generate a quadrature
@@ -726,7 +745,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[0]);
+    source_fe.generalized_face_support_points[face_no]);
 
   // Rule of thumb for FP accuracy,
   // that can be expected for a
@@ -745,11 +764,11 @@ FE_RaviartThomasNodal<dim>::get_subface_interpolation_matrix(
                                         0,
                                         subface);
 
-  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(); ++i)
+  for (unsigned int i = 0; i < source_fe.n_dofs_per_face(face_no); ++i)
     {
       const Point<dim> &p = subface_projection.point(i);
 
-      for (unsigned int j = 0; j < this->n_dofs_per_face(); ++j)
+      for (unsigned int j = 0; j < this->n_dofs_per_face(face_no); ++j)
         {
           double matrix_entry =
             this->shape_value_component(this->face_to_cell_index(j, 0), p, 0);
@@ -774,11 +793,11 @@ FE_RaviartThomasNodal<dim>::get_subface_interpolation_matrix(
   // this point. this must be so
   // since the shape functions sum up
   // to 1
-  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(); ++j)
+  for (unsigned int j = 0; j < source_fe.n_dofs_per_face(face_no); ++j)
     {
       double sum = 0.;
 
-      for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+      for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
         sum += interpolation_matrix(j, i);
 
       Assert(std::fabs(sum - 1) < 2e-13 * this->degree * (dim - 1),
diff --git a/source/fe/fe_rt_bubbles.cc b/source/fe/fe_rt_bubbles.cc
index 9f0b2025db..145ca3a659 100644
--- a/source/fe/fe_rt_bubbles.cc
+++ b/source/fe/fe_rt_bubbles.cc
@@ -75,18 +75,26 @@ FE_RT_Bubbles<dim>::FE_RT_Bubbles(const unsigned int deg)
       for (unsigned int i = 0; i < nc; ++i)
         this->prolongation[ref_case - 1][i].reinit(n_dofs, n_dofs);
     }
+
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   // Fill prolongation matrices with embedding operators
   // set tolerance to 1, as embedding error accumulate quickly
   FETools::compute_embedding_matrices(*this, this->prolongation, true, 1.0);
   FullMatrix<double> face_embeddings[GeometryInfo<dim>::max_children_per_face];
   for (unsigned int i = 0; i < GeometryInfo<dim>::max_children_per_face; ++i)
-    face_embeddings[i].reinit(this->n_dofs_per_face(), this->n_dofs_per_face());
+    face_embeddings[i].reinit(this->n_dofs_per_face(face_no),
+                              this->n_dofs_per_face(face_no));
   FETools::compute_face_embedding_matrices<dim, double>(*this,
                                                         face_embeddings,
                                                         0,
                                                         0);
-  this->interface_constraints.reinit((1 << (dim - 1)) * this->n_dofs_per_face(),
-                                     this->n_dofs_per_face());
+  this->interface_constraints.reinit((1 << (dim - 1)) *
+                                       this->n_dofs_per_face(face_no),
+                                     this->n_dofs_per_face(face_no));
   unsigned int target_row = 0;
   for (unsigned int d = 0; d < GeometryInfo<dim>::max_children_per_face; ++d)
     for (unsigned int i = 0; i < face_embeddings[d].m(); ++i)
@@ -131,8 +139,14 @@ template <int dim>
 void
 FE_RT_Bubbles<dim>::initialize_support_points(const unsigned int deg)
 {
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   this->generalized_support_points.resize(this->n_dofs_per_cell());
-  this->generalized_face_support_points[0].resize(this->n_dofs_per_face());
+  this->generalized_face_support_points[face_no].resize(
+    this->n_dofs_per_face(face_no));
 
   // Index of the point being entered
   unsigned int current = 0;
@@ -143,25 +157,28 @@ FE_RT_Bubbles<dim>::initialize_support_points(const unsigned int deg)
   if (dim > 1)
     {
       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[0][k] = face_points.point(k);
+      Assert(face_points.size() == this->n_dofs_per_face(face_no),
+             ExcInternalError());
+      for (unsigned int k = 0; k < this->n_dofs_per_face(face_no); ++k)
+        this->generalized_face_support_points[face_no][k] =
+          face_points.point(k);
       Quadrature<dim> faces =
         QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
                                               face_points);
-      for (unsigned int k = 0;
-           k < this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
+      for (unsigned int k = 0; k < this->n_dofs_per_face(face_no) *
+                                     GeometryInfo<dim>::faces_per_cell;
            ++k)
-        this->generalized_support_points[k] = faces.point(
-          k +
-          QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
-                                                   0,
-                                                   true,
-                                                   false,
-                                                   false,
-                                                   this->n_dofs_per_face()));
-
-      current = this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
+        this->generalized_support_points[k] =
+          faces.point(k + QProjector<dim>::DataSetDescriptor::face(
+                            this->reference_cell_type(),
+                            0,
+                            true,
+                            false,
+                            false,
+                            this->n_dofs_per_face(face_no)));
+
+      current =
+        this->n_dofs_per_face(face_no) * GeometryInfo<dim>::faces_per_cell;
     }
 
   if (deg == 1)
@@ -283,9 +300,9 @@ FE_RT_Bubbles<dim>::convert_generalized_support_point_values_to_dof_values(
   unsigned int fbase = 0;
   unsigned int f     = 0;
   for (; f < GeometryInfo<dim>::faces_per_cell;
-       ++f, fbase += this->n_dofs_per_face())
+       ++f, fbase += this->n_dofs_per_face(f))
     {
-      for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
+      for (unsigned int i = 0; i < this->n_dofs_per_face(f); ++i)
         {
           nodal_values[fbase + i] = support_point_values[fbase + i](
             GeometryInfo<dim>::unit_normal_direction[f]);
diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc
index 78d40b2ef0..de8169987e 100644
--- a/source/fe/fe_system.cc
+++ b/source/fe/fe_system.cc
@@ -858,7 +858,7 @@ FESystem<dim, spacedim>::face_to_cell_index(const unsigned int face_dof_index,
   // the DoFs within their own numbering. thus, translate to
   // the base element numbering and then back
   const std::pair<std::pair<unsigned int, unsigned int>, unsigned int>
-    face_base_index = this->face_system_to_base_index(face_dof_index);
+    face_base_index = this->face_system_to_base_index(face_dof_index, face);
 
   const unsigned int base_face_to_cell_index =
     this->base_element(face_base_index.first.first)
@@ -1392,6 +1392,11 @@ template <int dim, int spacedim>
 void
 FESystem<dim, spacedim>::build_interface_constraints()
 {
+  // TODO: the implementation makes the assumption that all faces have the
+  // same number of dofs
+  AssertDimension(this->n_unique_faces(), 1);
+  const unsigned int face_no = 0;
+
   // check whether all base elements implement their interface constraint
   // matrices. if this is not the case, then leave the interface costraints of
   // this composed element empty as well; however, the rest of the element is
@@ -1422,7 +1427,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[0][n];
+          n_index = this->face_system_to_base_table[face_no][n];
 
         // likewise for the m index. this is more complicated due to the
         // strange ordering we have for the dofs on the refined faces.
@@ -1461,7 +1466,7 @@ FESystem<dim, spacedim>::build_interface_constraints()
                     const unsigned int tmp1 =
                       2 * this->n_dofs_per_vertex() + index_in_line;
                     m_index.first =
-                      this->face_system_to_base_table[0][tmp1].first;
+                      this->face_system_to_base_table[face_no][tmp1].first;
 
                     // what we are still missing is the index of m within the
                     // base elements interface_constraints table
@@ -1472,12 +1477,13 @@ 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[0][tmp1].second >=
-                             2 * base_element(m_index.first.first)
-                                   .n_dofs_per_vertex(),
-                           ExcInternalError());
+                    Assert(
+                      this->face_system_to_base_table[face_no][tmp1].second >=
+                        2 *
+                          base_element(m_index.first.first).n_dofs_per_vertex(),
+                      ExcInternalError());
                     const unsigned int tmp2 =
-                      this->face_system_to_base_table[0][tmp1].second -
+                      this->face_system_to_base_table[face_no][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(),
@@ -1517,14 +1523,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[0][tmp1].first;
+                      this->face_system_to_base_table[face_no][tmp1].first;
 
-                    Assert(this->face_system_to_base_table[0][tmp1].second >=
-                             4 * base_element(m_index.first.first)
-                                   .n_dofs_per_vertex(),
-                           ExcInternalError());
+                    Assert(
+                      this->face_system_to_base_table[face_no][tmp1].second >=
+                        4 *
+                          base_element(m_index.first.first).n_dofs_per_vertex(),
+                      ExcInternalError());
                     const unsigned int tmp2 =
-                      this->face_system_to_base_table[0][tmp1].second -
+                      this->face_system_to_base_table[face_no][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(),
@@ -1543,42 +1550,45 @@ FESystem<dim, spacedim>::build_interface_constraints()
                     const unsigned int index_in_quad =
                       (m - 5 * this->n_dofs_per_vertex() -
                        12 * this->n_dofs_per_line()) %
-                      this->n_dofs_per_quad();
-                    Assert(index_in_quad < this->n_dofs_per_quad(),
+                      this->n_dofs_per_quad(face_no);
+                    Assert(index_in_quad < this->n_dofs_per_quad(face_no),
                            ExcInternalError());
                     const unsigned int sub_quad =
                       ((m - 5 * this->n_dofs_per_vertex() -
                         12 * this->n_dofs_per_line()) /
-                       this->n_dofs_per_quad());
+                       this->n_dofs_per_quad(face_no));
                     Assert(sub_quad < 4, ExcInternalError());
 
                     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[0].size(),
+                    Assert(tmp1 <
+                             this->face_system_to_base_table[face_no].size(),
                            ExcInternalError());
                     m_index.first =
-                      this->face_system_to_base_table[0][tmp1].first;
-
-                    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());
+                      this->face_system_to_base_table[face_no][tmp1].first;
+
+                    Assert(
+                      this->face_system_to_base_table[face_no][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[0][tmp1].second -
+                      this->face_system_to_base_table[face_no][tmp1].second -
                       4 *
                         base_element(m_index.first.first).n_dofs_per_vertex() -
                       4 * base_element(m_index.first.first).n_dofs_per_line();
                     Assert(tmp2 < base_element(m_index.first.first)
-                                    .n_dofs_per_quad(),
+                                    .n_dofs_per_quad(face_no),
                            ExcInternalError());
                     m_index.second =
                       5 *
                         base_element(m_index.first.first).n_dofs_per_vertex() +
                       12 * base_element(m_index.first.first).n_dofs_per_line() +
-                      base_element(m_index.first.first).n_dofs_per_quad() *
+                      base_element(m_index.first.first)
+                          .n_dofs_per_quad(face_no) *
                         sub_quad +
                       tmp2;
                   }
@@ -1648,18 +1658,25 @@ FESystem<dim, spacedim>::initialize(
     // FiniteElement
     this->system_to_component_table.resize(this->n_dofs_per_cell());
 
-    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,
                                             this->component_to_base_table,
                                             *this);
 
-    FETools::Compositing::build_face_tables(
-      this->face_system_to_base_table[0],
-      this->face_system_to_component_table[0],
-      *this);
+    this->face_system_to_component_table.resize(this->n_unique_faces());
+
+    for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no)
+      {
+        this->face_system_to_component_table[0].resize(
+          this->n_dofs_per_face(face_no));
+
+        FETools::Compositing::build_face_tables(
+          this->face_system_to_base_table[face_no],
+          this->face_system_to_component_table[face_no],
+          *this,
+          true,
+          face_no);
+      }
   }
 
   // now initialize interface constraints, support points, and other tables.
@@ -1702,43 +1719,56 @@ FESystem<dim, spacedim>::initialize(
   // initialize face support points (for dim==2,3). same procedure as above
   if (dim > 1)
     init_tasks += Threads::new_task([&]() {
-      // if one of the base elements has no support points, then it makes no
-      // sense to define support points for the composed element. In that case,
-      // return an empty array to demonstrate that fact (note that we ask
-      // whether the base element has no support points at all, not only none on
-      // the face!)
-      //
-      // on the other hand, if there is an element that simply has no degrees of
-      // freedom on the face at all, then we don't care whether it has support
-      // points or not. this is, for example, the case for the stable Stokes
-      // element Q(p)^dim \times DGP(p-1).
-      for (unsigned int base_el = 0; base_el < this->n_base_elements();
-           ++base_el)
-        if (!base_element(base_el).has_support_points() &&
-            (base_element(base_el).n_dofs_per_face() > 0))
-          {
-            this->unit_face_support_points[0].resize(0);
-            return;
-          }
-
+      for (unsigned int face_no = 0; face_no < this->n_unique_faces();
+           ++face_no)
+        {
+          // if one of the base elements has no support points, then it makes
+          // no sense to define support points for the composed element. In
+          // that case, return an empty array to demonstrate that fact (note
+          // that we ask whether the base element has no support points at
+          // all, not only none on the face!)
+          //
+          // on the other hand, if there is an element that simply has no
+          // degrees of freedom on the face at all, then we don't care whether
+          // it has support points or not. this is, for example, the case for
+          // the stable Stokes element Q(p)^dim \times DGP(p-1).
+          bool flag_has_no_support_points = false;
+
+          for (unsigned int base_el = 0; base_el < this->n_base_elements();
+               ++base_el)
+            if (!base_element(base_el).has_support_points() &&
+                (base_element(base_el).n_dofs_per_face(face_no) > 0))
+              {
+                this->unit_face_support_points[face_no].resize(0);
+                flag_has_no_support_points = true;
+                break;
+              }
 
-      // generate unit face support points from unit support points of sub
-      // elements
-      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[0][i].first.first;
-          const unsigned int index_in_base =
-            this->face_system_to_base_table[0][i].second;
+          if (flag_has_no_support_points)
+            continue;
 
-          Assert(index_in_base <
-                   base_element(base_i).unit_face_support_points[0].size(),
-                 ExcInternalError());
+          // generate unit face support points from unit support points of sub
+          // elements
+          this->unit_face_support_points[face_no].resize(
+            this->n_dofs_per_face(face_no));
 
-          this->unit_face_support_points[0][i] =
-            base_element(base_i).unit_face_support_points[0][index_in_base];
+          for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
+            {
+              const unsigned int base_i =
+                this->face_system_to_base_table[face_no][i].first.first;
+              const unsigned int index_in_base =
+                this->face_system_to_base_table[face_no][i].second;
+
+              Assert(
+                index_in_base <
+                  base_element(base_i).unit_face_support_points[face_no].size(),
+                ExcInternalError());
+
+              this->unit_face_support_points[face_no][i] =
+                base_element(base_i)
+                  .unit_face_support_points[face_no][index_in_base];
+            }
         }
     });
 
@@ -1813,36 +1843,40 @@ FESystem<dim, spacedim>::initialize(
   // initialize quad dof index permutation in 3d and higher
   if (dim >= 3)
     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[0]
-                 .n_elements() == 8 * this->n_dofs_per_quad(),
-             ExcInternalError());
-
-      // to obtain the shifts for this composed element, copy the shift
-      // information of the base elements
-      unsigned int index = 0;
-      for (unsigned int b = 0; b < this->n_base_elements(); ++b)
+      for (unsigned int face_no = 0; face_no < this->n_unique_faces();
+           ++face_no)
         {
-          const Table<2, int> &temp =
-            this->base_element(b)
-              .adjust_quad_dof_index_for_face_orientation_table[0];
-          for (unsigned int c = 0; c < this->element_multiplicity(b); ++c)
+          // the array into which we want to write should have the correct size
+          // already.
+          Assert(this->adjust_quad_dof_index_for_face_orientation_table[face_no]
+                     .n_elements() == 8 * this->n_dofs_per_quad(face_no),
+                 ExcInternalError());
+
+          // to obtain the shifts for this composed element, copy the shift
+          // information of the base elements
+          unsigned int index = 0;
+          for (unsigned int b = 0; b < this->n_base_elements(); ++b)
             {
-              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[0](
-                    index + i, j) = temp(i, j);
-              index += temp.size(0);
+              const Table<2, int> &temp =
+                this->base_element(b)
+                  .adjust_quad_dof_index_for_face_orientation_table[face_no];
+              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
+                        [face_no](index + i, j) = temp(i, j);
+                  index += temp.size(0);
+                }
             }
+          Assert(index == this->n_dofs_per_quad(face_no), ExcInternalError());
         }
-      Assert(index == this->n_dofs_per_quad(), ExcInternalError());
 
       // additionally compose the permutation information for lines
       Assert(this->adjust_line_dof_index_for_line_orientation_table.size() ==
                this->n_dofs_per_line(),
              ExcInternalError());
-      index = 0;
+      unsigned int index = 0;
       for (unsigned int b = 0; b < this->n_base_elements(); ++b)
         {
           const std::vector<int> &temp2 =
@@ -1887,12 +1921,12 @@ FESystem<dim, spacedim>::get_face_interpolation_matrix(
   FullMatrix<double> &                interpolation_matrix,
   const unsigned int                  face_no) const
 {
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+                              this->n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
 
   // since dofs for each base are independent, we only have to stack things up
   // from base element to base element
@@ -1928,8 +1962,8 @@ FESystem<dim, spacedim>::get_face_interpolation_matrix(
                  ExcNotImplemented());
 
           // get the interpolation from the bases
-          base_to_base_interpolation.reinit(base_other.n_dofs_per_face(),
-                                            base.n_dofs_per_face());
+          base_to_base_interpolation.reinit(base_other.n_dofs_per_face(face_no),
+                                            base.n_dofs_per_face(face_no));
           base.get_face_interpolation_matrix(base_other,
                                              base_to_base_interpolation,
                                              face_no);
@@ -1937,16 +1971,19 @@ FESystem<dim, spacedim>::get_face_interpolation_matrix(
           // now translate entries. we'd like to have something like
           // face_base_to_system_index, but that doesn't exist. rather, all we
           // have is the reverse. well, use that then
-          for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
-            if (this->face_system_to_base_index(i).first ==
+          for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
+            if (this->face_system_to_base_index(i, face_no).first ==
                 std::make_pair(base_index, multiplicity))
-              for (unsigned int j = 0; j < fe_other_system->n_dofs_per_face();
+              for (unsigned int j = 0;
+                   j < fe_other_system->n_dofs_per_face(face_no);
                    ++j)
-                if (fe_other_system->face_system_to_base_index(j).first ==
+                if (fe_other_system->face_system_to_base_index(j, face_no)
+                      .first ==
                     std::make_pair(base_index_other, multiplicity_other))
                   interpolation_matrix(j, i) = base_to_base_interpolation(
-                    fe_other_system->face_system_to_base_index(j).second,
-                    this->face_system_to_base_index(i).second);
+                    fe_other_system->face_system_to_base_index(j, face_no)
+                      .second,
+                    this->face_system_to_base_index(i, face_no).second);
 
           // advance to the next base element for this and the other fe_system;
           // see if we can simply advance the multiplicity by one, or if have to
@@ -2006,12 +2043,12 @@ FESystem<dim, spacedim>::get_subface_interpolation_matrix(
       (dynamic_cast<const FESystem<dim, spacedim> *>(&x_source_fe) != nullptr),
     (typename FiniteElement<dim, spacedim>::ExcInterpolationNotImplemented()));
 
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+                              this->n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
 
   // since dofs for each base are independent, we only have to stack things up
   // from base element to base element
@@ -2048,8 +2085,8 @@ FESystem<dim, spacedim>::get_subface_interpolation_matrix(
                  ExcNotImplemented());
 
           // get the interpolation from the bases
-          base_to_base_interpolation.reinit(base_other.n_dofs_per_face(),
-                                            base.n_dofs_per_face());
+          base_to_base_interpolation.reinit(base_other.n_dofs_per_face(face_no),
+                                            base.n_dofs_per_face(face_no));
           base.get_subface_interpolation_matrix(base_other,
                                                 subface,
                                                 base_to_base_interpolation,
@@ -2058,16 +2095,19 @@ FESystem<dim, spacedim>::get_subface_interpolation_matrix(
           // now translate entries. we'd like to have something like
           // face_base_to_system_index, but that doesn't exist. rather, all we
           // have is the reverse. well, use that then
-          for (unsigned int i = 0; i < this->n_dofs_per_face(); ++i)
-            if (this->face_system_to_base_index(i).first ==
+          for (unsigned int i = 0; i < this->n_dofs_per_face(face_no); ++i)
+            if (this->face_system_to_base_index(i, face_no).first ==
                 std::make_pair(base_index, multiplicity))
-              for (unsigned int j = 0; j < fe_other_system->n_dofs_per_face();
+              for (unsigned int j = 0;
+                   j < fe_other_system->n_dofs_per_face(face_no);
                    ++j)
-                if (fe_other_system->face_system_to_base_index(j).first ==
+                if (fe_other_system->face_system_to_base_index(j, face_no)
+                      .first ==
                     std::make_pair(base_index_other, multiplicity_other))
                   interpolation_matrix(j, i) = base_to_base_interpolation(
-                    fe_other_system->face_system_to_base_index(j).second,
-                    this->face_system_to_base_index(i).second);
+                    fe_other_system->face_system_to_base_index(j, face_no)
+                      .second,
+                    this->face_system_to_base_index(i, face_no).second);
 
           // advance to the next base element for this and the other fe_system;
           // see if we can simply advance the multiplicity by one, or if have to
@@ -2117,7 +2157,8 @@ template <int dim, int spacedim>
 template <int structdim>
 std::vector<std::pair<unsigned int, unsigned int>>
 FESystem<dim, spacedim>::hp_object_dof_identities(
-  const FiniteElement<dim, spacedim> &fe_other) const
+  const FiniteElement<dim, spacedim> &fe_other,
+  const unsigned int                  face_no) const
 {
   // since dofs on each subobject (vertex, line, ...) are ordered such that
   // first come all from the first base element all multiplicities, then
@@ -2167,7 +2208,8 @@ FESystem<dim, spacedim>::hp_object_dof_identities(
                 base_identities = base.hp_line_dof_identities(base_other);
                 break;
               case 2:
-                base_identities = base.hp_quad_dof_identities(base_other);
+                base_identities =
+                  base.hp_quad_dof_identities(base_other, face_no);
                 break;
               default:
                 Assert(false, ExcNotImplemented());
@@ -2247,9 +2289,9 @@ template <int dim, int spacedim>
 std::vector<std::pair<unsigned int, unsigned int>>
 FESystem<dim, spacedim>::hp_quad_dof_identities(
   const FiniteElement<dim, spacedim> &fe_other,
-  const unsigned int) const
+  const unsigned int                  face_no) const
 {
-  return hp_object_dof_identities<2>(fe_other);
+  return hp_object_dof_identities<2>(fe_other, face_no);
 }
 
 
@@ -2354,22 +2396,27 @@ FESystem<dim, spacedim>::unit_face_support_point(
   const unsigned int index,
   const unsigned int face_no) const
 {
-  AssertIndexRange(index, this->n_dofs_per_face());
-  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()));
+  AssertIndexRange(index, this->n_dofs_per_face(face_no));
+  Assert(
+    (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+       .size() == this->n_dofs_per_face(face_no)) ||
+      (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+         .size() == 0),
+    (typename FiniteElement<dim, spacedim>::ExcFEHasNoSupportPoints()));
 
   // let's see whether we have the information pre-computed
-  if (this->unit_face_support_points[0].size() != 0)
-    return this->unit_face_support_points[0][index];
+  if (this->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+        .size() != 0)
+    return this
+      ->unit_face_support_points[this->n_unique_faces() == 1 ? 0 : face_no]
+                                [index];
   else
     // no. ask the base element whether it would like to provide this
     // information
     return (
-      base_element(this->face_system_to_base_index(index).first.first)
-        .unit_face_support_point(this->face_system_to_base_index(index).second,
-                                 face_no));
+      base_element(this->face_system_to_base_index(index, face_no).first.first)
+        .unit_face_support_point(
+          this->face_system_to_base_index(index, face_no).second, face_no));
 }
 
 
diff --git a/source/fe/fe_tools.inst.in b/source/fe/fe_tools.inst.in
index c77d24dc53..981a5bc527 100644
--- a/source/fe/fe_tools.inst.in
+++ b/source/fe/fe_tools.inst.in
@@ -102,7 +102,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
           std::vector<std::pair<unsigned int, unsigned int>>
             &face_system_to_component_table,
           const FiniteElement<deal_II_dimension, deal_II_space_dimension> &fe,
-          const bool do_tensor_product);
+          const bool         do_tensor_product,
+          const unsigned int face_no);
 
       \}
 
diff --git a/source/fe/fe_trace.cc b/source/fe/fe_trace.cc
index 5190e0e655..584465e413 100644
--- a/source/fe/fe_trace.cc
+++ b/source/fe/fe_trace.cc
@@ -53,7 +53,10 @@ 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[0] = fe_q.get_unit_face_support_points();
+  AssertDimension(this->n_unique_faces(), fe_q.n_unique_faces());
+  for (unsigned int face_no = 0; face_no < this->n_unique_faces(); ++face_no)
+    this->unit_face_support_points[face_no] =
+      fe_q.get_unit_face_support_points(face_no);
 
   // 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
@@ -231,15 +234,15 @@ FE_TraceQ<dim, spacedim>::get_subface_interpolation_matrix(
   const FiniteElement<dim, spacedim> &x_source_fe,
   const unsigned int                  subface,
   FullMatrix<double> &                interpolation_matrix,
-  const unsigned int) const
+  const unsigned int                  face_no) const
 {
   // this is the code from FE_FaceQ
-  Assert(interpolation_matrix.n() == this->n_dofs_per_face(),
+  Assert(interpolation_matrix.n() == this->n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.n(),
-                              this->n_dofs_per_face()));
-  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(),
+                              this->n_dofs_per_face(face_no)));
+  Assert(interpolation_matrix.m() == x_source_fe.n_dofs_per_face(face_no),
          ExcDimensionMismatch(interpolation_matrix.m(),
-                              x_source_fe.n_dofs_per_face()));
+                              x_source_fe.n_dofs_per_face(face_no)));
 
   // see if source is a FaceQ element
   if (const FE_TraceQ<dim, spacedim> *source_fe =
@@ -247,7 +250,8 @@ FE_TraceQ<dim, spacedim>::get_subface_interpolation_matrix(
     {
       fe_q.get_subface_interpolation_matrix(source_fe->fe_q,
                                             subface,
-                                            interpolation_matrix);
+                                            interpolation_matrix,
+                                            face_no);
     }
   else if (dynamic_cast<const FE_Nothing<dim> *>(&x_source_fe) != nullptr)
     {