Make get_cell_range_category() more consistent with get_cell_category()

include/deal.II/matrix_free/matrix_free.h

@@ -1892,8 +1892,14 @@ class MatrixFree : public Subscriptor
   /**
    * Return the category the current batch range of cells was assigned to.
    * Categories run between the given values in the field
-   * AdditionalData::cell_vectorization_category for non-hp-DoFHandler types
+   * AdditionalData::cell_vectorization_category for the non-hp case
    * and return the active FE index in the hp-adaptive case.
+   *
+   * @note Following the behaviour of get_cell_category(), we return the
+   * maximum category of any cell batch. In the hp case, it is
+   * guaranteed that all cells and as a consequence all cell batches in a range
+   * have the same category. Otherwise, there may be different categories in
+   * different cell batches.
    */
   unsigned int
   get_cell_range_category(
@@ -1910,8 +1916,13 @@ class MatrixFree : public Subscriptor
   /**
    * Return the category the current batch of cells was assigned to. Categories
    * run between the given values in the field
-   * AdditionalData::cell_vectorization_category for non-hp-DoFHandler types
+   * AdditionalData::cell_vectorization_category for the non-hp case
    * and return the active FE index in the hp-adaptive case.
+   *
+   * @note In the non-hp case, a category of a cell batch is given
+   * as the maximum category of any of its cell. In the hp case or the case that
+   * MatrixFree::AdditionalData::cell_vectorization_categories_strict was
+   * enabled, it is guaranteed that all cells have the same category.
    */
   unsigned int
   get_cell_category(const unsigned int cell_batch_index) const;
@@ -2805,14 +2816,10 @@ inline unsigned int
 MatrixFree<dim, Number, VectorizedArrayType>::get_cell_range_category(
   const std::pair<unsigned int, unsigned int> range) const
 {
-  const auto result = get_cell_category(range.first);
+  auto result = get_cell_category(range.first);
 
-#  ifdef DEBUG
   for (unsigned int i = range.first; i < range.second; ++i)
-    Assert(result == get_cell_category(i),
-           ExcMessage(
-             "The cell batches of the range have different categories!"));
-#  endif
+    result = std::max(result, get_cell_category(i));
 
   return result;
 }
@@ -2824,14 +2831,13 @@ inline std::pair<unsigned int, unsigned int>
 MatrixFree<dim, Number, VectorizedArrayType>::get_face_range_category(
   const std::pair<unsigned int, unsigned int> range) const
 {
-  const auto result = get_face_category(range.first);
+  auto result = get_face_category(range.first);
 
-#  ifdef DEBUG
   for (unsigned int i = range.first; i < range.second; ++i)
-    Assert(result == get_face_category(i),
-           ExcMessage(
-             "The face batches of the range have different categories!"));
-#  endif
+    {
+      result.first  = std::max(result.first, get_face_category(i).first);
+      result.second = std::max(result.second, get_face_category(i).second);
+    }
 
   return result;
 }

tests/matrix_free/get_cell_catergory_01.cc

@@ -0,0 +1,246 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests MatrixFree::get_cell_category() and
+// MatrixFree::get_cell_range_category()
+
+#include <deal.II/distributed/shared_tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test()
+{
+  using Number              = double;
+  using VectorizedArrayType = VectorizedArray<Number>;
+
+  Triangulation<dim> tria;
+
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(3);
+
+  // caterorization - not strict
+  {
+    DoFHandler<dim> dof_handler(tria);
+    dof_handler.distribute_dofs(FE_Q<dim>(1));
+
+    typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, double>::AdditionalData::none;
+    data.mapping_update_flags  = update_quadrature_points;
+
+    std::vector<unsigned int> cell_vectorization_category(
+      tria.n_active_cells());
+    for (unsigned int i = 0; i < cell_vectorization_category.size(); ++i)
+      cell_vectorization_category[i] = i / 7;
+
+    data.cell_vectorization_category          = cell_vectorization_category;
+    data.cell_vectorization_categories_strict = false;
+
+    MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+    matrix_free.reinit(MappingQ1<dim>(),
+                       dof_handler,
+                       AffineConstraints<Number>(),
+                       QGauss<dim>(2),
+                       data);
+
+    double dummy = 0;
+
+    matrix_free.template cell_loop<double, double>(
+      [&](const auto &matrix_free, auto &, const auto &, const auto range) {
+        unsigned int max_category_range = 0;
+
+        for (unsigned int cell = range.first; cell < range.second; ++cell)
+          {
+            unsigned int max_category = 0;
+            for (unsigned int v = 0;
+                 v < matrix_free.n_active_entries_per_cell_batch(cell);
+                 ++v)
+              {
+                const auto category = cell_vectorization_category
+                  [matrix_free.get_cell_iterator(cell, v)->active_cell_index()];
+
+                // note: here and below we are using std::cout to output
+                // information that is useful for debugging; however, since the
+                // output depends on the hardware and might change if we
+                // internally modify the way we loop over cells, we don't write
+                // the data to the log but instead only check that the
+                // properties that are described in the documentation are
+                // fulfilled
+                std::cout << category << " ";
+                max_category = std::max(max_category, category);
+              }
+
+            std::cout << std::endl;
+
+            AssertDimension(max_category,
+                            matrix_free.get_cell_category(cell)); // test
+
+            max_category_range = std::max(max_category_range, max_category);
+          }
+
+        AssertDimension(max_category_range,
+                        matrix_free.get_cell_range_category(range)); // test
+        std::cout << std::endl;
+      },
+      dummy,
+      dummy);
+  }
+  deallog << "OK" << std::endl;
+
+  // caterorization - strict
+  {
+    DoFHandler<dim> dof_handler(tria);
+    dof_handler.distribute_dofs(FE_Q<dim>(1));
+
+    typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, double>::AdditionalData::none;
+    data.mapping_update_flags  = update_quadrature_points;
+    data.mapping_update_flags_boundary_faces = update_quadrature_points;
+    data.mapping_update_flags_inner_faces    = update_quadrature_points;
+    data.hold_all_faces_to_owned_cells       = true;
+
+    std::vector<unsigned int> cell_vectorization_category(
+      tria.n_active_cells());
+    for (unsigned int i = 0; i < cell_vectorization_category.size(); ++i)
+      cell_vectorization_category[i] = i / 7;
+
+    data.cell_vectorization_category          = cell_vectorization_category;
+    data.cell_vectorization_categories_strict = true;
+
+    MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+    matrix_free.reinit(MappingQ1<dim>(),
+                       dof_handler,
+                       AffineConstraints<Number>(),
+                       QGauss<dim>(2),
+                       data);
+
+    double dummy = 0;
+
+    matrix_free.template cell_loop<double, double>(
+      [&](const auto &matrix_free, auto &, const auto &, const auto range) {
+        unsigned int max_category_range = 0;
+
+        for (unsigned int cell = range.first; cell < range.second; ++cell)
+          {
+            for (unsigned int v = 0;
+                 v < matrix_free.n_active_entries_per_cell_batch(cell);
+                 ++v)
+              {
+                const auto category = cell_vectorization_category
+                  [matrix_free.get_cell_iterator(cell, v)->active_cell_index()];
+                std::cout << category << " ";
+                AssertDimension(category,
+                                matrix_free.get_cell_category(cell)); // test
+              }
+
+            std::cout << std::endl;
+
+            max_category_range =
+              std::max(max_category_range, matrix_free.get_cell_category(cell));
+          }
+
+        AssertDimension(max_category_range,
+                        matrix_free.get_cell_range_category(range)); // test
+        std::cout << std::endl;
+      },
+      dummy,
+      dummy);
+  }
+  deallog << "OK" << std::endl;
+
+
+  // hp
+  {
+    DoFHandler<dim> dof_handler(tria);
+
+    unsigned int counter = 0;
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        cell->set_active_fe_index((counter++) / 7);
+
+    hp::FECollection<dim> fe_collection;
+
+    for (unsigned int i = 0; i < (tria.n_cells() + 6) / 7; ++i)
+      fe_collection.push_back(FE_Q<dim>(1));
+
+    dof_handler.distribute_dofs(fe_collection);
+
+    typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, double>::AdditionalData::none;
+    data.mapping_update_flags  = update_quadrature_points;
+    data.mapping_update_flags_boundary_faces = update_quadrature_points;
+    data.mapping_update_flags_inner_faces    = update_quadrature_points;
+    data.hold_all_faces_to_owned_cells       = true;
+
+    MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+    matrix_free.reinit(MappingQ1<dim>(),
+                       dof_handler,
+                       AffineConstraints<Number>(),
+                       QGauss<dim>(2),
+                       data);
+
+
+    double dummy = 0;
+
+    matrix_free.template cell_loop<double, double>(
+      [&](const auto &matrix_free, auto &, const auto &, const auto range) {
+        unsigned int max_category_range = 0;
+
+        for (unsigned int cell = range.first; cell < range.second; ++cell)
+          {
+            for (unsigned int v = 0;
+                 v < matrix_free.n_active_entries_per_cell_batch(cell);
+                 ++v)
+              {
+                const auto category =
+                  matrix_free.get_cell_iterator(cell, v)->active_fe_index();
+                std::cout << category << " ";
+                AssertDimension(category,
+                                matrix_free.get_cell_category(cell)); // test
+                AssertDimension(
+                  category, matrix_free.get_cell_range_category(range)); // test
+              }
+
+            std::cout << std::endl;
+          }
+        std::cout << std::endl;
+      },
+      dummy,
+      dummy);
+  }
+  deallog << "OK" << std::endl;
+}
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  test<2>();
+}

tests/matrix_free/get_cell_catergory_01.output

@@ -0,0 +1,4 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK