/**
* 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(
/**
* 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;
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;
}
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;
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
projects / dealii.git / commitdiff