MatrixFree: use active_fe_indices within FaceComparator

diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h
index bf154269a244460b955dc30cab188a129e129f65..a2f0c5462e52ec3ad0d95f79114dbf71291858b3 100644 (file)
--- a/include/deal.II/matrix_free/face_setup_internal.h
+++ b/include/deal.II/matrix_free/face_setup_internal.h
@@ -1023,8 +1023,11 @@ namespace internal
      * to batch similar faces together for vectorization.
      */
     inline bool
-    compare_faces_for_vectorization(const FaceToCellTopology<1> &face1,
-                                    const FaceToCellTopology<1> &face2)
+    compare_faces_for_vectorization(
+      const FaceToCellTopology<1> &    face1,
+      const FaceToCellTopology<1> &    face2,
+      const std::vector<unsigned int> &active_fe_indices,
+      const unsigned int               length)
     {
       if (face1.interior_face_no != face2.interior_face_no)
         return false;
@@ -1034,6 +1037,19 @@ namespace internal
         return false;
       if (face1.face_orientation != face2.face_orientation)
         return false;
+
+      if (active_fe_indices.size() > 0)
+        {
+          if (active_fe_indices[face1.cells_interior[0] / length] !=
+              active_fe_indices[face2.cells_interior[0] / length])
+            return false;
+
+          if (face2.cells_exterior[0] != numbers::invalid_unsigned_int)
+            if (active_fe_indices[face1.cells_exterior[0] / length] !=
+                active_fe_indices[face2.cells_exterior[0] / length])
+              return false;
+        }
+
       return true;
     }
 
@@ -1048,10 +1064,37 @@ namespace internal
     template <int length>
     struct FaceComparator
     {
+      FaceComparator(const std::vector<unsigned int> &active_fe_indices)
+        : active_fe_indices(active_fe_indices)
+      {}
+
       bool
       operator()(const FaceToCellTopology<length> &face1,
                  const FaceToCellTopology<length> &face2) const
       {
+        // check if active fe indices match
+        if (active_fe_indices.size() > 0)
+          {
+            // ... for interior faces
+            if (active_fe_indices[face1.cells_interior[0] / length] <
+                active_fe_indices[face2.cells_interior[0] / length])
+              return true;
+            else if (active_fe_indices[face1.cells_interior[0] / length] >
+                     active_fe_indices[face2.cells_interior[0] / length])
+              return false;
+
+            // ... for exterior faces
+            if (face2.cells_exterior[0] != numbers::invalid_unsigned_int)
+              {
+                if (active_fe_indices[face1.cells_exterior[0] / length] <
+                    active_fe_indices[face2.cells_exterior[0] / length])
+                  return true;
+                else if (active_fe_indices[face1.cells_exterior[0] / length] >
+                         active_fe_indices[face2.cells_exterior[0] / length])
+                  return false;
+              }
+          }
+
         for (unsigned int i = 0; i < length; ++i)
           if (face1.cells_interior[i] < face2.cells_interior[i])
             return true;
@@ -1079,6 +1122,9 @@ namespace internal
 
         return false;
       }
+
+    private:
+      const std::vector<unsigned int> &active_fe_indices;
     };
 
 
@@ -1089,7 +1135,8 @@ namespace internal
       const std::vector<FaceToCellTopology<1>> &faces_in,
       const std::vector<bool> &                 hard_vectorization_boundary,
       std::vector<unsigned int> &               face_partition_data,
-      std::vector<FaceToCellTopology<vectorization_width>> &faces_out)
+      std::vector<FaceToCellTopology<vectorization_width>> &faces_out,
+      const std::vector<unsigned int> &                     active_fe_indices)
     {
       FaceToCellTopology<vectorization_width> macro_face;
       std::vector<std::vector<unsigned int>>  faces_type;
@@ -1119,7 +1166,9 @@ namespace internal
                 {
                   // Compare current face with first face of type type
                   if (compare_faces_for_vectorization(faces_in[face],
-                                                      faces_in[face_type[0]]))
+                                                      faces_in[face_type[0]],
+                                                      active_fe_indices,
+                                                      vectorization_width))
                     {
                       face_type.push_back(face);
                       goto face_found;
@@ -1131,9 +1180,11 @@ namespace internal
             }
 
           // insert new faces in sorted list to get good data locality
+          FaceComparator<vectorization_width> face_comparator(
+            active_fe_indices);
           std::set<FaceToCellTopology<vectorization_width>,
                    FaceComparator<vectorization_width>>
-            new_faces;
+            new_faces(face_comparator);
           for (const auto &face_type : faces_type)
             {
               macro_face.interior_face_no =

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 01809a8ea98ec3ec510a1fb003902ebbb9bcdf99..c0dce201685c589907385c14bc66089704f2a295 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -1403,7 +1403,8 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
         face_setup.inner_faces,
         hard_vectorization_boundary,
         task_info.face_partition_data,
-        face_info.faces);
+        face_info.faces,
+        dof_info[0].cell_active_fe_index);
 
       // on boundary faces, we must also respect the vectorization boundary of
       // the inner faces because we might have dependencies on ghosts of
@@ -1412,7 +1413,8 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
         face_setup.boundary_faces,
         hard_vectorization_boundary,
         task_info.boundary_partition_data,
-        face_info.faces);
+        face_info.faces,
+        dof_info[0].cell_active_fe_index);
 
       // for the other ghosted faces, there are no scheduling restrictions
       hard_vectorization_boundary.clear();
@@ -1422,7 +1424,8 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
         face_setup.inner_ghost_faces,
         hard_vectorization_boundary,
         task_info.ghost_face_partition_data,
-        face_info.faces);
+        face_info.faces,
+        dof_info[0].cell_active_fe_index);
       hard_vectorization_boundary.clear();
       hard_vectorization_boundary.resize(
         task_info.refinement_edge_face_partition_data.size(), false);
@@ -1430,7 +1433,8 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
         face_setup.refinement_edge_faces,
         hard_vectorization_boundary,
         task_info.refinement_edge_face_partition_data,
-        face_info.faces);
+        face_info.faces,
+        dof_info[0].cell_active_fe_index);
 
       cell_level_index.resize(
         cell_level_index.size() +