Fix MF::get_cell_active_fe_index() to recognise cell categorization

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 92a43ad4a856456519ea3eb9c982374255372ed9..81eef32ef60b3db2c91513e455255af303babb0a 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -7121,7 +7121,7 @@ inline FEEvaluation<dim,
                  dof_no,
                  quad_no,
                  first_selected_component,
-                 matrix_free.get_cell_active_fe_index(range))
+                 matrix_free.get_cell_active_fe_index(range, dof_no))
 {}
 
 

diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h
index e78ce79939caadba6963a1716e6667785bc89bcc..32cb3098d744eef0c4b7bac2f7738738a2ae4b75 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.h
+++ b/include/deal.II/matrix_free/matrix_free.h
@@ -1618,7 +1618,7 @@ public:
    */
   unsigned int
   get_cell_active_fe_index(
-    const std::pair<unsigned int, unsigned int> range) const;
+    const std::pair<unsigned int, unsigned int> range, unsigned int dof_index) const;
 
   /**
    * In the hp-adaptive case, return the active FE index of a face range.
@@ -2689,11 +2689,12 @@ MatrixFree<dim, Number, VectorizedArrayType>::n_active_fe_indices() const
 template <int dim, typename Number, typename VectorizedArrayType>
 unsigned int
 MatrixFree<dim, Number, VectorizedArrayType>::get_cell_active_fe_index(
-  const std::pair<unsigned int, unsigned int> range) const
+  const std::pair<unsigned int, unsigned int> range, unsigned int dof_handler_index) const
 {
   const auto &fe_indices = dof_info[0].cell_active_fe_index;
 
-  if (fe_indices.empty() == true)
+  if (fe_indices.empty() == true ||
+     dof_handlers[dof_handler_index]->get_fe_collection().size() == 1)
     return 0;
 
   const auto index = fe_indices[range.first];

diff --git a/tests/matrix_free/compute_diagonal_01.cc b/tests/matrix_free/compute_diagonal_01.cc
index 685838d8a72305c2c8028037edbaaf6be800ae4c..f9e2915f3c2bbba0810c3ac17c4e88978afef762 100644 (file)
--- a/tests/matrix_free/compute_diagonal_01.cc
+++ b/tests/matrix_free/compute_diagonal_01.cc
@@ -25,7 +25,7 @@ template <int dim,
           typename Number              = double,
           typename VectorizedArrayType = VectorizedArray<Number>>
 void
-test()
+test(bool use_categorization = false)
 {
   parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
 
@@ -38,6 +38,15 @@ test()
       cell->set_refine_flag();
   tria.execute_coarsening_and_refinement();
 
+  if (use_categorization)
+    {
+      for (auto &cell : tria.active_cell_iterators())
+        if (cell->is_active() && cell->is_locally_owned() &&
+            cell->center()[0] < 0.0)
+          cell->set_material_id(42);
+        else
+          cell->set_material_id(0);
+    }
   const FE_Q<dim>     fe_q(fe_degree);
   const FESystem<dim> fe(fe_q, n_components);
 
@@ -60,6 +69,18 @@ test()
     additional_data;
   additional_data.mapping_update_flags = update_values | update_gradients;
 
+  if (use_categorization)
+    {
+      additional_data.cell_vectorization_category.resize(tria.n_active_cells());
+      for (const auto &cell : tria.active_cell_iterators())
+        if (cell->is_locally_owned())
+          additional_data
+            .cell_vectorization_category[cell->active_cell_index()] =
+            cell->material_id();
+
+      additional_data.cell_vectorization_categories_strict = true;
+    }
+
   MappingQ<dim> mapping(1);
   QGauss<1>     quad(fe_degree + 1);
 
@@ -96,4 +117,9 @@ main(int argc, char **argv)
 
   test<2, 1, 2, 1, float>(); // scalar
   test<2, 1, 2, 2, float>(); // vector
+
+  // Same as above, but testing additional categorization of vectorized cell
+  // batches
+  test<2, 1, 2, 1, float>(true); // scalar
+  test<2, 1, 2, 2, float>(true); // vector
 }