From a3b7cea558105a1d139b606573468e2089d36ea2 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Thu, 4 Mar 2021 18:02:53 +0100
Subject: [PATCH] MatrixFree: introduce FaceToCellTopology::face_type

---
 .../deal.II/matrix_free/evaluation_kernels.h  |   4 +-
 include/deal.II/matrix_free/face_info.h       |   6 +
 .../deal.II/matrix_free/face_setup_internal.h |  15 ++
 include/deal.II/matrix_free/fe_evaluation.h   |  74 +++++---
 .../matrix_free/mapping_info.templates.h      |   2 -
 .../matrix_free/matrix_free.templates.h       | 171 ++++++++++++------
 .../matrix_free/shape_info.templates.h        |   7 +-
 7 files changed, 189 insertions(+), 90 deletions(-)

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 12abbb4ccd..a9ea830a82 100644
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -2549,7 +2549,7 @@ namespace internal
       if (data.element_type == MatrixFreeFunctions::tensor_none)
         {
           const unsigned int n_dofs     = data.dofs_per_component_on_cell;
-          const unsigned int n_q_points = data.n_q_points_face;
+          const unsigned int n_q_points = data.n_q_points_faces[face_no];
           const auto         shape_info = data.data.front();
 
           using Eval = EvaluatorTensorProduct<evaluate_general,
@@ -2704,7 +2704,7 @@ namespace internal
       if (data.element_type == MatrixFreeFunctions::tensor_none)
         {
           const unsigned int n_dofs     = data.dofs_per_component_on_cell;
-          const unsigned int n_q_points = data.n_q_points_face;
+          const unsigned int n_q_points = data.n_q_points_faces[face_no];
           const auto         shape_info = data.data.front();
 
           using Eval = EvaluatorTensorProduct<evaluate_general,
diff --git a/include/deal.II/matrix_free/face_info.h b/include/deal.II/matrix_free/face_info.h
index b411c27d7e..2e1e01ef5a 100644
--- a/include/deal.II/matrix_free/face_info.h
+++ b/include/deal.II/matrix_free/face_info.h
@@ -110,6 +110,12 @@ namespace internal
        */
       unsigned char face_orientation;
 
+      /**
+       * Reference-cell type of the given face: 0 for line or quadrilateral,
+       * 1 for triangle.
+       */
+      unsigned char face_type;
+
       /**
        * Return the memory consumption of the present data structure.
        */
diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h
index d480a4db25..e3905e6fad 100644
--- a/include/deal.II/matrix_free/face_setup_internal.h
+++ b/include/deal.II/matrix_free/face_setup_internal.h
@@ -780,6 +780,9 @@ namespace internal
                         info.cells_exterior[0] = numbers::invalid_unsigned_int;
                         info.interior_face_no  = f;
                         info.exterior_face_no  = dcell->face(f)->boundary_id();
+                        info.face_type =
+                          dcell->face(f)->reference_cell() !=
+                          dealii::ReferenceCells::get_hypercube<dim - 1>();
                         info.subface_index =
                           GeometryInfo<dim>::max_children_per_cell;
                         info.face_orientation = 0;
@@ -959,6 +962,9 @@ namespace internal
       else
         info.exterior_face_no = cell->neighbor_face_no(face_no);
 
+      info.face_type = cell->face(face_no)->reference_cell() !=
+                       dealii::ReferenceCells::get_hypercube<dim - 1>();
+
       info.subface_index = GeometryInfo<dim>::max_children_per_cell;
       Assert(neighbor->level() <= cell->level(), ExcInternalError());
       if (cell->level() > neighbor->level())
@@ -1037,6 +1043,8 @@ namespace internal
         return false;
       if (face1.face_orientation != face2.face_orientation)
         return false;
+      if (face1.face_type != face2.face_type)
+        return false;
 
       if (active_fe_indices.size() > 0)
         {
@@ -1072,6 +1080,12 @@ namespace internal
       operator()(const FaceToCellTopology<length> &face1,
                  const FaceToCellTopology<length> &face2) const
       {
+        // check if active FE indices match
+        if (face1.face_type < face2.face_type)
+          return true;
+        else if (face1.face_type > face2.face_type)
+          return false;
+
         // check if active FE indices match
         if (active_fe_indices.size() > 0)
           {
@@ -1187,6 +1201,7 @@ namespace internal
             new_faces(face_comparator);
           for (const auto &face_type : faces_type)
             {
+              macro_face.face_type = faces_in[face_type[0]].face_type;
               macro_face.interior_face_no =
                 faces_in[face_type[0]].interior_face_no;
               macro_face.exterior_face_no =
diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 45de596e71..1dd0f581d6 100644
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -298,7 +298,8 @@ protected:
     const unsigned int n_q_points,
     const bool         is_interior_face,
     const unsigned int active_fe_index,
-    const unsigned int active_quad_index);
+    const unsigned int active_quad_index,
+    const unsigned int face_type);
 
   /**
    * Constructor that comes with reduced functionality and works similar as
@@ -1143,7 +1144,8 @@ protected:
     const unsigned int n_q_points,
     const bool         is_interior_face,
     const unsigned int active_fe_index,
-    const unsigned int active_quad_index);
+    const unsigned int active_quad_index,
+    const unsigned int face_type);
 
   /**
    * Constructor that comes with reduced functionality and works similar as
@@ -1428,7 +1430,8 @@ protected:
     const unsigned int n_q_points,
     const bool         is_interior_face  = true,
     const unsigned int active_fe_index   = numbers::invalid_unsigned_int,
-    const unsigned int active_quad_index = numbers::invalid_unsigned_int);
+    const unsigned int active_quad_index = numbers::invalid_unsigned_int,
+    const unsigned int face_type         = numbers::invalid_unsigned_int);
 
   /**
    * Constructor with reduced functionality for similar usage of FEEvaluation
@@ -1565,7 +1568,8 @@ protected:
     const unsigned int n_q_points,
     const bool         is_interior_face  = true,
     const unsigned int active_fe_index   = numbers::invalid_unsigned_int,
-    const unsigned int active_quad_index = numbers::invalid_unsigned_int);
+    const unsigned int active_quad_index = numbers::invalid_unsigned_int,
+    const unsigned int face_type         = numbers::invalid_unsigned_int);
 
   /**
    * Constructor with reduced functionality for similar usage of FEEvaluation
@@ -1726,7 +1730,8 @@ protected:
     const unsigned int n_q_points,
     const bool         is_interior_face  = true,
     const unsigned int active_fe_index   = numbers::invalid_unsigned_int,
-    const unsigned int active_quad_index = numbers::invalid_unsigned_int);
+    const unsigned int active_quad_index = numbers::invalid_unsigned_int,
+    const unsigned int face_type         = numbers::invalid_unsigned_int);
 
   /**
    * Constructor with reduced functionality for similar usage of FEEvaluation
@@ -1872,7 +1877,8 @@ protected:
     const unsigned int                                n_q_points,
     const bool                                        is_interior_face = true,
     const unsigned int active_fe_index   = numbers::invalid_unsigned_int,
-    const unsigned int active_quad_index = numbers::invalid_unsigned_int);
+    const unsigned int active_quad_index = numbers::invalid_unsigned_int,
+    const unsigned int face_type         = numbers::invalid_unsigned_int);
 
   /**
    * Constructor with reduced functionality for similar usage of FEEvaluation
@@ -3069,6 +3075,9 @@ public:
    * DoFHandler/AffineConstraints pair the given evaluator should be attached
    * to.
    *
+   * @param face_type In the case of a face, indicate its reference-cell type
+   * (0 for line or quadrilateral 1 for triangle).
+   *
    * @param active_quad_index If matrix_free was set up with hp::Collection
    * objects, this parameter selects the appropriate number of the quadrature
    * formula.
@@ -3080,7 +3089,8 @@ public:
     const unsigned int                                  quad_no          = 0,
     const unsigned int first_selected_component                          = 0,
     const unsigned int active_fe_index   = numbers::invalid_unsigned_int,
-    const unsigned int active_quad_index = numbers::invalid_unsigned_int);
+    const unsigned int active_quad_index = numbers::invalid_unsigned_int,
+    const unsigned int face_type         = numbers::invalid_unsigned_int);
 
   /**
    * Constructor. Takes all data stored in MatrixFree for a given face range,
@@ -3342,7 +3352,8 @@ inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
     const unsigned int n_q_points,
     const bool         is_interior_face,
     const unsigned int active_fe_index_in,
-    const unsigned int active_quad_index_in)
+    const unsigned int active_quad_index_in,
+    const unsigned int face_type)
   : scratch_data_array(data_in.acquire_scratch_data())
   , quad_no(quad_no_in)
   , matrix_info(&data_in)
@@ -3366,7 +3377,12 @@ inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
            active_quad_index_in :
            std::min<unsigned int>(active_fe_index,
                                   mapping_data->descriptor.size() - 1)))
-  , descriptor(&mapping_data->descriptor[active_quad_index])
+  , descriptor(
+      &mapping_data->descriptor
+         [is_face ?
+            (active_quad_index * std::max<unsigned int>(1, dim - 1) +
+             (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
+            active_quad_index])
   , n_quadrature_points(descriptor->n_q_points)
   , data(&data_in.get_shape_info(
       dof_no,
@@ -4019,7 +4035,8 @@ inline FEEvaluationBase<dim,
                    const unsigned int n_q_points,
                    const bool         is_interior_face,
                    const unsigned int active_fe_index,
-                   const unsigned int active_quad_index)
+                   const unsigned int active_quad_index,
+                   const unsigned int face_type)
   : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
       data_in,
       dof_no,
@@ -4029,7 +4046,8 @@ inline FEEvaluationBase<dim,
       n_q_points,
       is_interior_face,
       active_fe_index,
-      active_quad_index)
+      active_quad_index,
+      face_type)
   , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
   , dof_values_initialized(false)
   , values_quad_initialized(false)
@@ -6159,7 +6177,8 @@ inline FEEvaluationAccess<dim,
     const unsigned int n_q_points,
     const bool         is_interior_face,
     const unsigned int active_fe_index,
-    const unsigned int active_quad_index)
+    const unsigned int active_quad_index,
+    const unsigned int face_type)
   : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
       data_in,
       dof_no,
@@ -6169,7 +6188,8 @@ inline FEEvaluationAccess<dim,
       n_q_points,
       is_interior_face,
       active_fe_index,
-      active_quad_index)
+      active_quad_index,
+      face_type)
 {}
 
 
@@ -6265,7 +6285,8 @@ inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
     const unsigned int n_q_points,
     const bool         is_interior_face,
     const unsigned int active_fe_index,
-    const unsigned int active_quad_index)
+    const unsigned int active_quad_index,
+    const unsigned int face_type)
   : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
       data_in,
       dof_no,
@@ -6275,7 +6296,8 @@ inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
       n_q_points,
       is_interior_face,
       active_fe_index,
-      active_quad_index)
+      active_quad_index,
+      face_type)
 {}
 
 
@@ -6576,7 +6598,8 @@ inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
     const unsigned int n_q_points,
     const bool         is_interior_face,
     const unsigned int active_fe_index,
-    const unsigned int active_quad_index)
+    const unsigned int active_quad_index,
+    const unsigned int face_type)
   : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
       data_in,
       dof_no,
@@ -6586,7 +6609,8 @@ inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
       n_q_points,
       is_interior_face,
       active_fe_index,
-      active_quad_index)
+      active_quad_index,
+      face_type)
 {}
 
 
@@ -6985,7 +7009,8 @@ inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
                      const unsigned int n_q_points,
                      const bool         is_interior_face,
                      const unsigned int active_fe_index,
-                     const unsigned int active_quad_index)
+                     const unsigned int active_quad_index,
+                     const unsigned int face_type)
   : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
       data_in,
       dof_no,
@@ -6995,7 +7020,8 @@ inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
       n_q_points,
       is_interior_face,
       active_fe_index,
-      active_quad_index)
+      active_quad_index,
+      face_type)
 {}
 
 
@@ -8412,7 +8438,8 @@ inline FEFaceEvaluation<dim,
     const unsigned int                                  quad_no,
     const unsigned int first_selected_component,
     const unsigned int active_fe_index,
-    const unsigned int active_quad_index)
+    const unsigned int active_quad_index,
+    const unsigned int face_type)
   : BaseClass(matrix_free,
               dof_no,
               first_selected_component,
@@ -8421,7 +8448,8 @@ inline FEFaceEvaluation<dim,
               static_n_q_points,
               is_interior_face,
               active_fe_index,
-              active_quad_index)
+              active_quad_index,
+              face_type)
   , dofs_per_component(this->data->dofs_per_component_on_cell)
   , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
   , n_q_points(this->n_quadrature_points)
@@ -8454,7 +8482,9 @@ inline FEFaceEvaluation<dim,
                      quad_no,
                      first_selected_component,
                      matrix_free.get_face_active_fe_index(range,
-                                                          is_interior_face))
+                                                          is_interior_face),
+                     numbers::invalid_unsigned_int,
+                     matrix_free.get_face_info(range.first).face_type)
 {}
 
 
diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index d51a855ba7..445e28dea7 100644
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -1843,8 +1843,6 @@ namespace internal
 #ifndef DEAL_II_WITH_SIMPLEX_SUPPORT
               // currently only non-hp-case...
               AssertDimension(mapping_in.size(), 1);
-              AssertDimension(mapping_info.face_data[my_q].descriptor.size(),
-                              1);
 #endif
 
               // We assume that we have the faces sorted by the active FE
diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 4229b31ba3..147b71f572 100644
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -134,33 +134,45 @@ namespace
   {
   public:
     FaceRangeCompartor(const std::vector<unsigned int> &fe_indices,
-                       const bool                       include)
+                       const bool                       include,
+                       const bool                       only_face_type)
       : fe_indices(fe_indices)
       , include(include)
+      , only_face_type(only_face_type)
     {}
 
     template <int vectorization_width>
     bool
     operator()(const internal::MatrixFreeFunctions::FaceToCellTopology<
-                 vectorization_width> &face,
-               const unsigned int &    fe_index)
+                 vectorization_width> &           face,
+               const std::array<unsigned int, 2> &fe_index)
     {
-      const unsigned int fe_index_face =
-        fe_indices[face.cells_interior[0] / vectorization_width];
+      const unsigned int face_type = face.face_type;
 
-      return fe_index_face < fe_index ||
-             (include ? (fe_index_face == fe_index) : false);
+      if (only_face_type)
+        return include ? (face_type <= fe_index[0]) : (face_type < fe_index[0]);
+
+      const std::array<unsigned int, 2> fe_index_face = {
+        {face_type, fe_indices[face.cells_interior[0] / vectorization_width]}};
+
+      return include ? (fe_index_face <= fe_index) : (fe_index_face < fe_index);
     }
 
     template <int vectorization_width>
     bool
     operator()(const internal::MatrixFreeFunctions::FaceToCellTopology<
-                 vectorization_width> &                     face,
-               const std::pair<unsigned int, unsigned int> &fe_index)
+                 vectorization_width> &           face,
+               const std::array<unsigned int, 3> &fe_index)
     {
-      const std::pair<unsigned int, unsigned int> fe_index_face = {
-        fe_indices[face.cells_interior[0] / vectorization_width],
-        fe_indices[face.cells_exterior[0] / vectorization_width]};
+      const unsigned int face_type = face.face_type;
+
+      if (only_face_type)
+        return include ? (face_type <= fe_index[0]) : (face_type < fe_index[0]);
+
+      const std::array<unsigned int, 3> fe_index_face = {
+        {face_type,
+         fe_indices[face.cells_interior[0] / vectorization_width],
+         fe_indices[face.cells_exterior[0] / vectorization_width]}};
 
       return include ? (fe_index_face <= fe_index) : (fe_index_face < fe_index);
     }
@@ -168,6 +180,7 @@ namespace
   private:
     const std::vector<unsigned int> &fe_indices;
     const bool                       include;
+    const bool                       only_face_type;
   };
 } // namespace
 
@@ -540,13 +553,18 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
 
           const auto create_inner_face_subrange_hp_by_index =
             [&](const std::pair<unsigned int, unsigned int> &range,
+                const unsigned int                           face_type,
                 const unsigned int                           fe_index_interior,
                 const unsigned int                           fe_index_exterior,
                 const unsigned int                           dof_handler_index =
                   0) -> std::pair<unsigned int, unsigned int> {
-            if (dof_info[dof_handler_index].max_fe_index == 0 ||
-                dof_handlers[dof_handler_index]->get_fe_collection().size() ==
-                  1)
+            const unsigned int n_face_types =
+              std::max<unsigned int>(dim - 1, 1);
+
+            if ((dof_info[dof_handler_index].max_fe_index == 0 ||
+                 dof_handlers[dof_handler_index]->get_fe_collection().size() ==
+                   1) &&
+                n_face_types == 1)
               return range;
 
             AssertIndexRange(fe_index_interior,
@@ -555,24 +573,32 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
                              dof_info[dof_handler_index].max_fe_index);
             const std::vector<unsigned int> &fe_indices =
               dof_info[dof_handler_index].cell_active_fe_index;
-            if (fe_indices.empty() == true)
+            if (fe_indices.empty() == true && n_face_types == 1)
               return range;
             else
               {
+                const bool only_face_type =
+                  dof_info[dof_handler_index].max_fe_index == 0 ||
+                  dof_handlers[dof_handler_index]->get_fe_collection().size() ==
+                    1 ||
+                  fe_indices.empty() == true;
+
                 std::pair<unsigned int, unsigned int> return_range;
                 return_range.first =
-                  std::lower_bound(face_info.faces.begin() + range.first,
-                                   face_info.faces.begin() + range.second,
-                                   std::pair<unsigned int, unsigned int>{
-                                     fe_index_interior, fe_index_exterior},
-                                   FaceRangeCompartor(fe_indices, false)) -
+                  std::lower_bound(
+                    face_info.faces.begin() + range.first,
+                    face_info.faces.begin() + range.second,
+                    std::array<unsigned int, 3>{
+                      {face_type, fe_index_interior, fe_index_exterior}},
+                    FaceRangeCompartor(fe_indices, false, only_face_type)) -
                   face_info.faces.begin();
                 return_range.second =
-                  std::lower_bound(face_info.faces.begin() + return_range.first,
-                                   face_info.faces.begin() + range.second,
-                                   std::pair<unsigned int, unsigned int>{
-                                     fe_index_interior, fe_index_exterior},
-                                   FaceRangeCompartor(fe_indices, true)) -
+                  std::lower_bound(
+                    face_info.faces.begin() + return_range.first,
+                    face_info.faces.begin() + range.second,
+                    std::array<unsigned int, 3>{
+                      {face_type, fe_index_interior, fe_index_exterior}},
+                    FaceRangeCompartor(fe_indices, true, only_face_type)) -
                   face_info.faces.begin();
                 Assert(return_range.first >= range.first &&
                          return_range.second <= range.second,
@@ -592,20 +618,27 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
                     task_info.face_partition_data[i + 1]};
 
                   if (range.second > range.first)
-                    for (unsigned int i = 0; i < this->n_active_fe_indices();
-                         ++i)
-                      for (unsigned int j = 0; j < this->n_active_fe_indices();
-                           ++j)
-                        {
-                          const auto subrange =
-                            create_inner_face_subrange_hp_by_index(range, i, j);
+                    for (unsigned int t = 0;
+                         t < std::max<unsigned int>(dim - 1, 1);
+                         ++t)
+                      for (unsigned int i = 0; i < this->n_active_fe_indices();
+                           ++i)
+                        for (unsigned int j = 0;
+                             j < this->n_active_fe_indices();
+                             ++j)
+                          {
+                            const auto subrange =
+                              create_inner_face_subrange_hp_by_index(range,
+                                                                     t,
+                                                                     i,
+                                                                     j);
 
-                          if (subrange.second <= subrange.first)
-                            continue;
+                            if (subrange.second <= subrange.first)
+                              continue;
 
-                          data.push_back(subrange.first);
-                          data.push_back(subrange.second);
-                        }
+                            data.push_back(subrange.first);
+                            data.push_back(subrange.second);
+                          }
                 }
 
               ptr.push_back(data.size() / 2);
@@ -624,34 +657,47 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
 
           const auto create_boundary_face_subrange_hp_by_index =
             [&](const std::pair<unsigned int, unsigned int> &range,
+                const unsigned int                           face_type,
                 const unsigned int                           fe_index,
                 const unsigned int                           dof_handler_index =
                   0) -> std::pair<unsigned int, unsigned int> {
-            if (dof_info[dof_handler_index].max_fe_index == 0 ||
-                dof_handlers[dof_handler_index]->get_fe_collection().size() ==
-                  1)
+            const unsigned int n_face_types =
+              std::max<unsigned int>(dim - 1, 1);
+
+            if ((dof_info[dof_handler_index].max_fe_index == 0 ||
+                 dof_handlers[dof_handler_index]->get_fe_collection().size() ==
+                   1) &&
+                n_face_types == 1)
               return range;
 
             AssertIndexRange(fe_index,
                              dof_info[dof_handler_index].max_fe_index);
             const std::vector<unsigned int> &fe_indices =
               dof_info[dof_handler_index].cell_active_fe_index;
-            if (fe_indices.empty() == true)
+            if (fe_indices.empty() == true && n_face_types == 1)
               return range;
             else
               {
+                const bool only_face_type =
+                  dof_info[dof_handler_index].max_fe_index == 0 ||
+                  dof_handlers[dof_handler_index]->get_fe_collection().size() ==
+                    1 ||
+                  fe_indices.empty() == true;
+
                 std::pair<unsigned int, unsigned int> return_range;
                 return_range.first =
-                  std::lower_bound(face_info.faces.begin() + range.first,
-                                   face_info.faces.begin() + range.second,
-                                   fe_index,
-                                   FaceRangeCompartor(fe_indices, false)) -
+                  std::lower_bound(
+                    face_info.faces.begin() + range.first,
+                    face_info.faces.begin() + range.second,
+                    std::array<unsigned int, 2>{{face_type, fe_index}},
+                    FaceRangeCompartor(fe_indices, false, only_face_type)) -
                   face_info.faces.begin();
                 return_range.second =
-                  std::lower_bound(face_info.faces.begin() + return_range.first,
-                                   face_info.faces.begin() + range.second,
-                                   fe_index,
-                                   FaceRangeCompartor(fe_indices, true)) -
+                  std::lower_bound(
+                    face_info.faces.begin() + return_range.first,
+                    face_info.faces.begin() + range.second,
+                    std::array<unsigned int, 2>{{face_type, fe_index}},
+                    FaceRangeCompartor(fe_indices, true, only_face_type)) -
                   face_info.faces.begin();
                 Assert(return_range.first >= range.first &&
                          return_range.second <= range.second,
@@ -672,18 +718,23 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
                     task_info.boundary_partition_data[i + 1]};
 
                   if (range.second > range.first)
-                    for (unsigned int i = 0; i < this->n_active_fe_indices();
-                         ++i)
-                      {
-                        const auto cell_subrange =
-                          create_boundary_face_subrange_hp_by_index(range, i);
+                    for (unsigned int t = 0;
+                         t < std::max<unsigned int>(dim - 1, 1);
+                         ++t)
+                      for (unsigned int i = 0; i < this->n_active_fe_indices();
+                           ++i)
+                        {
+                          const auto cell_subrange =
+                            create_boundary_face_subrange_hp_by_index(range,
+                                                                      t,
+                                                                      i);
 
-                        if (cell_subrange.second <= cell_subrange.first)
-                          continue;
+                          if (cell_subrange.second <= cell_subrange.first)
+                            continue;
 
-                        data.push_back(cell_subrange.first);
-                        data.push_back(cell_subrange.second);
-                      }
+                          data.push_back(cell_subrange.first);
+                          data.push_back(cell_subrange.second);
+                        }
                 }
 
               ptr.push_back(data.size() / 2);
diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index c3293a0f12..c489a86041 100644
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -93,10 +93,9 @@ namespace internal
 
       const auto fe_poly = dynamic_cast<const FE_Poly<dim, dim> *>(&fe);
 
-      if (dynamic_cast<const FE_SimplexP<dim, dim> *>(&fe) != nullptr ||
-          dynamic_cast<const FE_SimplexDGP<dim, dim> *>(&fe) != nullptr ||
-          dynamic_cast<const FE_WedgeP<dim, dim> *>(&fe) != nullptr ||
-          dynamic_cast<const FE_PyramidP<dim, dim> *>(&fe) != nullptr)
+      if (dynamic_cast<const FE_SimplexPoly<dim, dim> *>(&fe) != nullptr ||
+          dynamic_cast<const FE_Wedge<dim, dim> *>(&fe) != nullptr ||
+          dynamic_cast<const FE_Pyramid<dim, dim> *>(&fe) != nullptr)
         {
           scalar_lexicographic.resize(fe.n_dofs_per_cell());
           for (unsigned int i = 0; i < scalar_lexicographic.size(); ++i)
-- 
2.39.5