Make MappingInfo::initialize_faces work for multiple fe indices

diff --git a/include/deal.II/matrix_free/mapping_info.h b/include/deal.II/matrix_free/mapping_info.h
index 941d530556351b7809cfcf01eba7f483081223d0..dc63926076e7b8c136188270aad631bfcade4488 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info.h
+++ b/include/deal.II/matrix_free/mapping_info.h
@@ -503,6 +503,7 @@ namespace internal
         const std::vector<std::pair<unsigned int, unsigned int>> &cells,
         const std::vector<FaceToCellTopology<VectorizedArrayType::size()>>
           &                                       faces,
+        const std::vector<unsigned int> &         active_fe_index,
         const dealii::hp::MappingCollection<dim> &mapping);
 
       /**

diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h
index 2a669655d560917d854e52ecc3ba5008dad8ca8e..6bea2049f1393c7774441e7e2212465a2cbde612 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info.templates.h
+++ b/include/deal.II/matrix_free/mapping_info.templates.h
@@ -450,7 +450,8 @@ namespace internal
           // Could call these functions in parallel, but not useful because
           // the work inside is nicely split up already
           initialize_cells(tria, cells, active_fe_index, *mapping);
-          initialize_faces(tria, cells, face_info.faces, *mapping);
+          initialize_faces(
+            tria, cells, face_info.faces, active_fe_index, *mapping);
           initialize_faces_by_cells(tria, cells, *mapping);
         }
     }
@@ -487,7 +488,8 @@ namespace internal
           // Could call these functions in parallel, but not useful because
           // the work inside is nicely split up already
           initialize_cells(tria, cells, active_fe_index, *mapping);
-          initialize_faces(tria, cells, face_info.faces, *mapping);
+          initialize_faces(
+            tria, cells, face_info.faces, active_fe_index, *mapping);
           initialize_faces_by_cells(tria, cells, *mapping);
         }
     }
@@ -889,14 +891,13 @@ namespace internal
         std::vector<std::vector<std::shared_ptr<dealii::FEValues<dim>>>>
           fe_values(mapping_info.cell_data.size());
 
-        const unsigned int n_active_fe_indices =
+        const unsigned int max_active_fe_index =
           active_fe_index.size() > 0 ?
-            (*std::max_element(active_fe_index.begin(), active_fe_index.end()) +
-             1) :
-            1;
+            *std::max_element(active_fe_index.begin(), active_fe_index.end()) :
+            0;
 
         for (unsigned int i = 0; i < fe_values.size(); ++i)
-          fe_values[i].resize(n_active_fe_indices);
+          fe_values[i].resize(max_active_fe_index + 1);
         const UpdateFlags update_flags = mapping_info.update_flags_cells;
         const UpdateFlags update_flags_feval =
           (update_flags & update_jacobians ? update_jacobians :
@@ -1759,6 +1760,7 @@ namespace internal
         const std::vector<std::pair<unsigned int, unsigned int>> &cells,
         const std::vector<FaceToCellTopology<VectorizedArrayType::size()>>
           &                                            faces,
+        const std::vector<unsigned int> &              active_fe_index,
         const dealii::hp::MappingCollection<dim> &     mapping_in,
         MappingInfo<dim, Number, VectorizedArrayType> &mapping_info,
         std::pair<
@@ -1766,24 +1768,25 @@ namespace internal
             MappingInfoStorage<dim - 1, dim, Number, VectorizedArrayType>>,
           CompressedFaceData<dim, Number, VectorizedArrayType>> &data)
       {
-        AssertDimension(mapping_in.size(), 1);
-        const auto &mapping = mapping_in[0];
-
         FE_Nothing<dim> dummy_fe;
 
+        const unsigned int max_active_fe_index =
+          active_fe_index.size() > 0 ?
+            *std::max_element(active_fe_index.begin(), active_fe_index.end()) :
+            0;
+
         std::vector<std::vector<std::shared_ptr<FEFaceValues<dim>>>>
           fe_face_values_container(mapping_info.face_data.size());
         for (unsigned int my_q = 0; my_q < mapping_info.face_data.size();
              ++my_q)
-          fe_face_values_container[my_q].resize(
-            mapping_info.face_data[my_q].descriptor.size());
+          fe_face_values_container[my_q].resize(max_active_fe_index + 1);
 
         std::vector<std::vector<std::shared_ptr<FEFaceValues<dim>>>>
           fe_boundary_face_values_container(mapping_info.face_data.size());
         for (unsigned int my_q = 0; my_q < mapping_info.face_data.size();
              ++my_q)
-          fe_boundary_face_values_container[my_q].resize(
-            mapping_info.face_data[my_q].descriptor.size());
+          fe_boundary_face_values_container[my_q].resize(max_active_fe_index +
+                                                         1);
 
         std::vector<std::vector<std::shared_ptr<FESubfaceValues<dim>>>>
           fe_subface_values_container(mapping_info.face_data.size());
@@ -1801,25 +1804,47 @@ namespace internal
           for (unsigned int my_q = 0; my_q < mapping_info.face_data.size();
                ++my_q)
             {
+#ifndef DEAL_II_WITH_SIMPLEX_SUPPORT
               // currently only non-hp case...
-              Assert(mapping_info.face_data[my_q].descriptor.size() == 1,
-                     ExcNotImplemented());
+              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
+              // indices so that the active fe index of the interior side of the
+              // face batch is the same as the fe index of the interior side of
+              // its first entry.
+              const unsigned int fe_index =
+                active_fe_index.size() > 0 ?
+                  active_fe_index[faces[face].cells_interior[0] /
+                                  VectorizedArrayType::size()] :
+                  0;
+              const unsigned int hp_quad_index =
+                mapping_info.cell_data[my_q].descriptor.size() == 1 ? 0 :
+                                                                      fe_index;
+              const unsigned int hp_mapping_index =
+                mapping_in.size() == 1 ? 0 : fe_index;
+
+              const auto &               mapping = mapping_in[hp_mapping_index];
               const Quadrature<dim - 1> &quadrature =
-                mapping_info.face_data[my_q].descriptor[0].quadrature;
+                mapping_info.face_data[my_q]
+                  .descriptor[hp_quad_index]
+                  .quadrature;
 
               const bool is_boundary_face =
                 faces[face].cells_exterior[0] == numbers::invalid_unsigned_int;
 
               if (is_boundary_face &&
-                  fe_boundary_face_values_container[my_q][0] == nullptr)
-                fe_boundary_face_values_container[my_q][0] =
+                  fe_boundary_face_values_container[my_q][fe_index] == nullptr)
+                fe_boundary_face_values_container[my_q][fe_index] =
                   std::make_shared<FEFaceValues<dim>>(
                     mapping,
                     dummy_fe,
                     quadrature,
                     mapping_info.update_flags_boundary_faces);
-              else if (fe_face_values_container[my_q][0] == nullptr)
-                fe_face_values_container[my_q][0] =
+              else if (fe_face_values_container[my_q][fe_index] == nullptr)
+                fe_face_values_container[my_q][fe_index] =
                   std::make_shared<FEFaceValues<dim>>(
                     mapping,
                     dummy_fe,
@@ -1827,8 +1852,9 @@ namespace internal
                     mapping_info.update_flags_inner_faces);
 
               FEFaceValues<dim> &fe_face_values =
-                is_boundary_face ? *fe_boundary_face_values_container[my_q][0] :
-                                   *fe_face_values_container[my_q][0];
+                is_boundary_face ?
+                  *fe_boundary_face_values_container[my_q][fe_index] :
+                  *fe_face_values_container[my_q][fe_index];
               const unsigned int n_q_points =
                 fe_face_values.n_quadrature_points;
               face_data.resize(n_q_points);
@@ -1962,9 +1988,44 @@ namespace internal
                       if (faces[face].subface_index >=
                           GeometryInfo<dim>::max_children_per_cell)
                         {
-                          fe_face_values.reinit(cell_it,
-                                                faces[face].exterior_face_no);
-                          actual_fe_face_values = &fe_face_values;
+                          // We assume that we have the faces sorted by the
+                          // active fe indices so that the active fe index of
+                          // the exterior side of the face batch is the same as
+                          // the fe index of the exterior side of its first
+                          // entry.
+                          const unsigned int fe_index =
+                            active_fe_index.size() > 0 ?
+                              active_fe_index[faces[face].cells_exterior[0] /
+                                              VectorizedArrayType::size()] :
+                              0;
+                          const unsigned int hp_quad_index =
+                            mapping_info.cell_data[my_q].descriptor.size() ==
+                                1 ?
+                              0 :
+                              fe_index;
+                          const unsigned int hp_mapping_index =
+                            mapping_in.size() == 1 ? 0 : fe_index;
+
+                          const auto &mapping = mapping_in[hp_mapping_index];
+                          const Quadrature<dim - 1> &quadrature =
+                            mapping_info.face_data[my_q]
+                              .descriptor[hp_quad_index]
+                              .quadrature;
+
+                          if (fe_face_values_container[my_q][fe_index] ==
+                              nullptr)
+                            fe_face_values_container[my_q][fe_index] =
+                              std::make_shared<FEFaceValues<dim>>(
+                                mapping,
+                                dummy_fe,
+                                quadrature,
+                                mapping_info.update_flags_boundary_faces);
+
+                          fe_face_values_container[my_q][fe_index]->reinit(
+                            cell_it, faces[face].exterior_face_no);
+
+                          actual_fe_face_values =
+                            fe_face_values_container[my_q][fe_index].get();
                         }
                       else
                         {
@@ -2406,6 +2467,7 @@ namespace internal
       const dealii::Triangulation<dim> &                                  tria,
       const std::vector<std::pair<unsigned int, unsigned int>> &          cells,
       const std::vector<FaceToCellTopology<VectorizedArrayType::size()>> &faces,
+      const std::vector<unsigned int> &         active_fe_index,
       const dealii::hp::MappingCollection<dim> &mapping)
     {
       face_type.resize(faces.size(), general);
@@ -2449,6 +2511,7 @@ namespace internal
               tria,
               cells,
               faces,
+              active_fe_index,
               mapping,
               *this,
               data_faces_local.back());