]> https://gitweb.dealii.org/ - dealii.git/commitdiff
MatrixFree: split up ranges according to active_fe_indices 11260/head
authorPeter Munch <peterrmuench@gmail.com>
Tue, 1 Dec 2020 18:54:26 +0000 (19:54 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Wed, 2 Dec 2020 14:22:43 +0000 (15:22 +0100)
include/deal.II/matrix_free/matrix_free.h
include/deal.II/matrix_free/matrix_free.templates.h
tests/simplex/matrix_free_02.cc
tests/simplex/matrix_free_range_iteration_01.cc [new file with mode: 0644]
tests/simplex/matrix_free_range_iteration_01.output [new file with mode: 0644]

index 4438f81211e556e693ff90667e4f6c05d9c379cb..e04a41b761316de6a43510550ee17f1449678e6c 100644 (file)
@@ -1481,6 +1481,51 @@ public:
     const unsigned int                           fe_index,
     const unsigned int                           dof_handler_index = 0) const;
 
+  /**
+   * In the hp adaptive case, a subrange of internal faces as computed during
+   * loop() might contain internal faces with elements of different active
+   * fe indices. Use this function to compute what the subrange for a given pair
+   * of active fe indices is.
+   */
+  std::pair<unsigned int, unsigned int>
+  create_inner_face_subrange_hp_by_index(
+    const std::pair<unsigned int, unsigned int> &range,
+    const unsigned int                           fe_index_interior,
+    const unsigned int                           fe_index_exterior,
+    const unsigned int                           dof_handler_index = 0) const;
+
+  /**
+   * In the hp adaptive case, a subrange of boundary faces as computed during
+   * loop() might contain boundary faces with elements of different active
+   * fe indices. Use this function to compute what the subrange for a given
+   * active fe indices is.
+   */
+  std::pair<unsigned int, unsigned int>
+  create_boundary_face_subrange_hp_by_index(
+    const std::pair<unsigned int, unsigned int> &range,
+    const unsigned int                           fe_index,
+    const unsigned int                           dof_handler_index = 0) const;
+
+  /**
+   * In the hp adaptive case, return number of active_fe_indices.
+   */
+  unsigned int
+  n_active_fe_indices() const;
+
+  /**
+   * In the hp adaptive case, return the active_fe_index of a cell range.
+   */
+  unsigned int
+  get_cell_active_fe_index(
+    const std::pair<unsigned int, unsigned int> range) const;
+
+  /**
+   * In the hp adaptive case, return the active_fe_index of a face range.
+   */
+  unsigned int
+  get_face_active_fe_index(const std::pair<unsigned int, unsigned int> range,
+                           const bool is_interior_face = true) const;
+
   //@}
 
   /**
@@ -2471,6 +2516,75 @@ MatrixFree<dim, Number, VectorizedArrayType>::at_irregular_cell(
 
 
 
+template <int dim, typename Number, typename VectorizedArrayType>
+unsigned int
+MatrixFree<dim, Number, VectorizedArrayType>::n_active_fe_indices() const
+{
+  return shape_info.size(2);
+}
+
+
+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 auto &fe_indices = dof_info[0].cell_active_fe_index;
+
+  if (fe_indices.empty() == true)
+    return 0;
+
+  const auto index = fe_indices[range.first];
+
+  for (unsigned int i = range.first; i < range.second; ++i)
+    AssertDimension(index, fe_indices[i]);
+
+  return index;
+}
+
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+unsigned int
+MatrixFree<dim, Number, VectorizedArrayType>::get_face_active_fe_index(
+  const std::pair<unsigned int, unsigned int> range,
+  const bool                                  is_interior_face) const
+{
+  const auto &fe_indices = dof_info[0].cell_active_fe_index;
+
+  if (fe_indices.empty() == true)
+    return 0;
+
+  if (is_interior_face)
+    {
+      const unsigned int index =
+        fe_indices[face_info.faces[range.first].cells_interior[0] /
+                   VectorizedArrayType::size()];
+
+      for (unsigned int i = range.first; i < range.second; ++i)
+        AssertDimension(index,
+                        fe_indices[face_info.faces[i].cells_interior[0] /
+                                   VectorizedArrayType::size()]);
+
+      return index;
+    }
+  else
+    {
+      const unsigned int index =
+        fe_indices[face_info.faces[range.first].cells_exterior[0] /
+                   VectorizedArrayType::size()];
+
+      for (unsigned int i = range.first; i < range.second; ++i)
+        AssertDimension(index,
+                        fe_indices[face_info.faces[i].cells_exterior[0] /
+                                   VectorizedArrayType::size()]);
+
+      return index;
+    }
+}
+
+
+
 template <int dim, typename Number, typename VectorizedArrayType>
 inline unsigned int
 MatrixFree<dim, Number, VectorizedArrayType>::n_components_filled(
@@ -4525,8 +4639,17 @@ namespace internal
     cell(const std::pair<unsigned int, unsigned int> &cell_range) override
     {
       if (cell_function != nullptr && cell_range.second > cell_range.first)
-        (container.*
-         cell_function)(matrix_free, this->dst, this->src, cell_range);
+        for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
+          {
+            const auto cell_subrange =
+              matrix_free.create_cell_subrange_hp_by_index(cell_range, i);
+
+            if (cell_subrange.second <= cell_subrange.first)
+              continue;
+
+            (container.*
+             cell_function)(matrix_free, this->dst, this->src, cell_subrange);
+          }
     }
 
     // Runs the assembler on interior faces. If no function is given, nothing
@@ -4535,8 +4658,19 @@ namespace internal
     face(const std::pair<unsigned int, unsigned int> &face_range) override
     {
       if (face_function != nullptr && face_range.second > face_range.first)
-        (container.*
-         face_function)(matrix_free, this->dst, this->src, face_range);
+        for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
+          for (unsigned int j = 0; j < matrix_free.n_active_fe_indices(); ++j)
+            {
+              const auto face_subrange =
+                matrix_free.create_inner_face_subrange_hp_by_index(face_range,
+                                                                   i,
+                                                                   j);
+
+              if (face_subrange.second <= face_subrange.first)
+                continue;
+              (container.*
+               face_function)(matrix_free, this->dst, this->src, face_subrange);
+            }
     }
 
     // Runs the assembler on boundary faces. If no function is given, nothing
@@ -4545,8 +4679,20 @@ namespace internal
     boundary(const std::pair<unsigned int, unsigned int> &face_range) override
     {
       if (boundary_function != nullptr && face_range.second > face_range.first)
-        (container.*
-         boundary_function)(matrix_free, this->dst, this->src, face_range);
+        for (unsigned int i = 0; i < matrix_free.n_active_fe_indices(); ++i)
+          {
+            const auto face_subrange =
+              matrix_free.create_boundary_face_subrange_hp_by_index(face_range,
+                                                                    i);
+
+            if (face_subrange.second <= face_subrange.first)
+              continue;
+
+            (container.*boundary_function)(matrix_free,
+                                           this->dst,
+                                           this->src,
+                                           face_subrange);
+          }
     }
 
     // Starts the communication for the update ghost values operation. We
index e2c0e006f36f28de026ac90090503e6e740ea239..de3b9ef39a1d11bc52dfebc717cf3e584d1a1b50 100644 (file)
@@ -85,12 +85,17 @@ std::pair<unsigned int, unsigned int>
 MatrixFree<dim, Number, VectorizedArrayType>::create_cell_subrange_hp_by_index(
   const std::pair<unsigned int, unsigned int> &range,
   const unsigned int                           fe_index,
-  const unsigned int                           vector_component) const
+  const unsigned int                           dof_handler_index) const
 {
-  AssertIndexRange(fe_index, dof_info[vector_component].max_fe_index);
+  if (dof_info[dof_handler_index].max_fe_index == 0)
+    return range;
+
+  AssertIndexRange(fe_index, dof_info[dof_handler_index].max_fe_index);
   const std::vector<unsigned int> &fe_indices =
-    dof_info[vector_component].cell_active_fe_index;
-  if (fe_indices.empty() == true)
+    dof_info[dof_handler_index].cell_active_fe_index;
+
+  if (fe_indices.empty() == true ||
+      dof_handlers[dof_handler_index]->get_fe_collection().size() == 1)
     return range;
   else
     {
@@ -124,14 +129,145 @@ MatrixFree<dim, Number, VectorizedArrayType>::create_cell_subrange_hp_by_index(
 
 
 
+namespace
+{
+  class FaceRangeCompartor
+  {
+  public:
+    FaceRangeCompartor(const std::vector<unsigned int> &fe_indices,
+                       const bool                       include)
+      : fe_indices(fe_indices)
+      , include(include)
+    {}
+
+    template <int vectorization_width>
+    bool
+    operator()(const internal::MatrixFreeFunctions::FaceToCellTopology<
+                 vectorization_width> &face,
+               const unsigned int &    fe_index)
+    {
+      const unsigned int fe_index_face =
+        fe_indices[face.cells_interior[0] / vectorization_width];
+
+      return fe_index_face < fe_index ||
+             (include ? (fe_index_face == fe_index) : false);
+    }
+
+    template <int vectorization_width>
+    bool
+    operator()(const internal::MatrixFreeFunctions::FaceToCellTopology<
+                 vectorization_width> &                     face,
+               const std::pair<unsigned int, unsigned int> &fe_index)
+    {
+      const std::pair<unsigned int, unsigned int> fe_index_face = {
+        fe_indices[face.cells_interior[0] / vectorization_width],
+        fe_indices[face.cells_exterior[0] / vectorization_width]};
+
+      return include ? (fe_index_face <= fe_index) : (fe_index_face < fe_index);
+    }
+
+  private:
+    const std::vector<unsigned int> &fe_indices;
+    const bool                       include;
+  };
+} // namespace
+
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+std::pair<unsigned int, unsigned int>
+MatrixFree<dim, Number, VectorizedArrayType>::
+  create_inner_face_subrange_hp_by_index(
+    const std::pair<unsigned int, unsigned int> &range,
+    const unsigned int                           fe_index_interior,
+    const unsigned int                           fe_index_exterior,
+    const unsigned int                           dof_handler_index) const
+{
+  if (dof_info[dof_handler_index].max_fe_index == 0 ||
+      dof_handlers[dof_handler_index]->get_fe_collection().size() == 1)
+    return range;
+
+  AssertIndexRange(fe_index_interior, dof_info[dof_handler_index].max_fe_index);
+  AssertIndexRange(fe_index_exterior, dof_info[dof_handler_index].max_fe_index);
+  const std::vector<unsigned int> &fe_indices =
+    dof_info[dof_handler_index].cell_active_fe_index;
+  if (fe_indices.empty() == true)
+    return range;
+  else
+    {
+      std::pair<unsigned int, unsigned int> return_range;
+      return_range.first =
+        std::lower_bound(face_info.faces.begin() + range.first,
+                         face_info.faces.begin() + range.second,
+                         std::pair<unsigned int, unsigned int>{
+                           fe_index_interior, fe_index_exterior},
+                         FaceRangeCompartor(fe_indices, false)) -
+        face_info.faces.begin();
+      return_range.second =
+        std::lower_bound(face_info.faces.begin() + return_range.first,
+                         face_info.faces.begin() + range.second,
+                         std::pair<unsigned int, unsigned int>{
+                           fe_index_interior, fe_index_exterior},
+                         FaceRangeCompartor(fe_indices, true)) -
+        face_info.faces.begin();
+      Assert(return_range.first >= range.first &&
+               return_range.second <= range.second,
+             ExcInternalError());
+      return return_range;
+    }
+}
+
+
+
+template <int dim, typename Number, typename VectorizedArrayType>
+std::pair<unsigned int, unsigned int>
+MatrixFree<dim, Number, VectorizedArrayType>::
+  create_boundary_face_subrange_hp_by_index(
+    const std::pair<unsigned int, unsigned int> &range,
+    const unsigned int                           fe_index,
+    const unsigned int                           dof_handler_index) const
+{
+  if (dof_info[dof_handler_index].max_fe_index == 0 ||
+      dof_handlers[dof_handler_index]->get_fe_collection().size() == 1)
+    return range;
+
+  AssertIndexRange(fe_index, dof_info[dof_handler_index].max_fe_index);
+  const std::vector<unsigned int> &fe_indices =
+    dof_info[dof_handler_index].cell_active_fe_index;
+  if (fe_indices.empty() == true)
+    return range;
+  else
+    {
+      std::pair<unsigned int, unsigned int> return_range;
+      return_range.first =
+        std::lower_bound(face_info.faces.begin() + range.first,
+                         face_info.faces.begin() + range.second,
+                         fe_index,
+                         FaceRangeCompartor(fe_indices, false)) -
+        face_info.faces.begin();
+      return_range.second =
+        std::lower_bound(face_info.faces.begin() + return_range.first,
+                         face_info.faces.begin() + range.second,
+                         fe_index,
+                         FaceRangeCompartor(fe_indices, true)) -
+        face_info.faces.begin();
+      Assert(return_range.first >= range.first &&
+               return_range.second <= range.second,
+             ExcInternalError());
+      return return_range;
+    }
+}
+
+
+
 template <int dim, typename Number, typename VectorizedArrayType>
 void
 MatrixFree<dim, Number, VectorizedArrayType>::renumber_dofs(
   std::vector<types::global_dof_index> &renumbering,
-  const unsigned int                    vector_component)
+  const unsigned int                    dof_handler_index)
 {
-  AssertIndexRange(vector_component, dof_info.size());
-  dof_info[vector_component].compute_dof_renumbering(renumbering);
+  AssertIndexRange(dof_handler_index, dof_info.size());
+  dof_info[dof_handler_index].compute_dof_renumbering(renumbering);
 }
 
 
@@ -1546,7 +1682,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
 #ifdef DEAL_II_WITH_MPI
   {
     // non-buffering mode is only supported if the indices of all cells are
-    // contiguous for all dof_info objects and hp is not enabled.
+    // contiguous for all dof_info objects.
     bool is_non_buffering_sm_supported = true;
     for (const auto &di : dof_info)
       {
index 7cb9ade094a2692c4bb7a501731ad148e36f31ab..d884c35976f0301bf62c7b63511f02ad9b623a06 100644 (file)
 
 #include "../tests.h"
 
-#include "./simplex_grids.h"
-
 using namespace dealii;
 
+namespace dealii
+{
+  namespace GridGenerator
+  {
+    template <int dim, int spacedim>
+    void
+    subdivided_hyper_rectangle_with_simplices_mix(
+      Triangulation<dim, spacedim> &   tria,
+      const std::vector<unsigned int> &repetitions,
+      const Point<dim> &               p1,
+      const Point<dim> &               p2,
+      const bool                       colorize = false)
+    {
+      AssertDimension(dim, spacedim);
+
+      AssertThrow(colorize == false, ExcNotImplemented());
+
+      std::vector<Point<spacedim>> vertices;
+      std::vector<CellData<dim>>   cells;
+
+      if (dim == 2)
+        {
+          // determine cell sizes
+          const Point<dim> dx((p2[0] - p1[0]) / repetitions[0],
+                              (p2[1] - p1[1]) / repetitions[1]);
+
+          // create vertices
+          for (unsigned int j = 0; j <= repetitions[1]; ++j)
+            for (unsigned int i = 0; i <= repetitions[0]; ++i)
+              vertices.push_back(
+                Point<spacedim>(p1[0] + dx[0] * i, p1[1] + dx[1] * j));
+
+          // create cells
+          for (unsigned int j = 0; j < repetitions[1]; ++j)
+            for (unsigned int i = 0; i < repetitions[0]; ++i)
+              {
+                // create reference QUAD cell
+                std::array<unsigned int, 4> quad{{
+                  (j + 0) * (repetitions[0] + 1) + i + 0, //
+                  (j + 0) * (repetitions[0] + 1) + i + 1, //
+                  (j + 1) * (repetitions[0] + 1) + i + 0, //
+                  (j + 1) * (repetitions[0] + 1) + i + 1  //
+                }};                                       //
+
+                if (j < repetitions[1] / 2 && i < repetitions[0] / 2)
+                  {
+                    CellData<dim> quad_;
+                    quad_.vertices = {quad[0], quad[1], quad[2], quad[3]};
+                    cells.push_back(quad_);
+
+                    continue;
+                  }
+
+                // TRI cell 0
+                {
+                  CellData<dim> tri;
+                  tri.vertices = {quad[0], quad[1], quad[2]};
+                  cells.push_back(tri);
+                }
+
+                // TRI cell 1
+                {
+                  CellData<dim> tri;
+                  tri.vertices = {quad[3], quad[2], quad[1]};
+                  cells.push_back(tri);
+                }
+              }
+        }
+      else
+        {
+          AssertThrow(colorize == false, ExcNotImplemented());
+        }
+
+      // actually create triangulation
+      tria.create_triangulation(vertices, cells, SubCellData());
+    }
+
+
+    template <int dim, int spacedim>
+    void
+    subdivided_hyper_cube_with_simplices_mix(Triangulation<dim, spacedim> &tria,
+                                             const unsigned int repetitions,
+                                             const double       p1 = 0.0,
+                                             const double       p2 = 1.0,
+                                             const bool colorize   = false)
+    {
+      if (dim == 2)
+        {
+          subdivided_hyper_rectangle_with_simplices_mix(
+            tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2}, colorize);
+        }
+      else if (dim == 3)
+        {
+          subdivided_hyper_rectangle_with_simplices_mix(
+            tria,
+            {{repetitions, repetitions, repetitions}},
+            {p1, p1, p1},
+            {p2, p2, p2},
+            colorize);
+        }
+      else
+        {
+          AssertThrow(false, ExcNotImplemented())
+        }
+    }
+  } // namespace GridGenerator
+} // namespace dealii
+
 template <int dim>
 class PoissonOperator
 {
 public:
-  using VectorType = LinearAlgebra::distributed::Vector<double>;
+  using VectorType       = LinearAlgebra::distributed::Vector<double>;
+  using FECellIntegrator = FEEvaluation<dim, -1, 0, 1, double>;
 
   PoissonOperator(const MatrixFree<dim, double> &matrix_free,
                   const bool                     do_helmholtz)
@@ -75,24 +182,17 @@ public:
     const int dummy = 0;
 
     matrix_free.template cell_loop<VectorType, int>(
-      [&](const auto &data, auto &dst, const auto &, const auto cell_range) {
-        for (unsigned int i = 0; i < 2; ++i)
-          {
-            const auto cell_subrange =
-              data.create_cell_subrange_hp_by_index(cell_range, i);
+      [&](const auto &data, auto &dst, const auto &, const auto range) {
+        const auto       i = data.get_cell_active_fe_index(range);
+        FECellIntegrator phi(matrix_free, 0, 0, 0, i, i);
 
-            FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free, 0, 0, 0, i, i);
-
-            for (unsigned int cell = cell_subrange.first;
-                 cell < cell_subrange.second;
-                 ++cell)
-              {
-                phi.reinit(cell);
-                for (unsigned int q = 0; q < phi.n_q_points; ++q)
-                  phi.submit_value(1.0, q);
+        for (unsigned int cell = range.first; cell < range.second; ++cell)
+          {
+            phi.reinit(cell);
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
+              phi.submit_value(1.0, q);
 
-                phi.integrate_scatter(true, false, dst);
-              }
+            phi.integrate_scatter(true, false, dst);
           }
       },
       vec,
@@ -105,30 +205,24 @@ public:
   vmult(VectorType &dst, const VectorType &src) const
   {
     matrix_free.template cell_loop<VectorType, VectorType>(
-      [&](const auto &data, auto &dst, const auto &src, const auto cell_range) {
-        for (unsigned int i = 0; i < 2; ++i)
+      [&](const auto &data, auto &dst, const auto &src, const auto range) {
+        const auto       i = data.get_cell_active_fe_index(range);
+        FECellIntegrator phi(matrix_free, 0, 0, 0, i, i);
+
+        for (unsigned int cell = range.first; cell < range.second; ++cell)
           {
-            const auto cell_subrange =
-              data.create_cell_subrange_hp_by_index(cell_range, i);
+            phi.reinit(cell);
+            phi.gather_evaluate(src, do_helmholtz, true);
 
-            FEEvaluation<dim, -1, 0, 1, double> phi(matrix_free, 0, 0, 0, i, i);
-            for (unsigned int cell = cell_subrange.first;
-                 cell < cell_subrange.second;
-                 ++cell)
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
               {
-                phi.reinit(cell);
-                phi.gather_evaluate(src, do_helmholtz, true);
+                if (do_helmholtz)
+                  phi.submit_value(phi.get_value(q), q);
 
-                for (unsigned int q = 0; q < phi.n_q_points; ++q)
-                  {
-                    if (do_helmholtz)
-                      phi.submit_value(phi.get_value(q), q);
-
-                    phi.submit_gradient(phi.get_gradient(q), q);
-                  }
-
-                phi.integrate_scatter(do_helmholtz, true, dst);
+                phi.submit_gradient(phi.get_gradient(q), q);
               }
+
+            phi.integrate_scatter(do_helmholtz, true, dst);
           }
       },
       dst,
diff --git a/tests/simplex/matrix_free_range_iteration_01.cc b/tests/simplex/matrix_free_range_iteration_01.cc
new file mode 100644 (file)
index 0000000..54479d8
--- /dev/null
@@ -0,0 +1,164 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test ShapeData for Simplex::FE_P and Simplex::QGauss
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <ios>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test()
+{
+  const unsigned int degree = 1;
+
+  Triangulation<dim> tria;
+  GridGenerator::subdivided_hyper_cube(tria, 4);
+
+  hp::FECollection<dim> fe{FE_Q<2>(degree), FE_Q<2>(degree)};
+  MappingQGeneric<dim>  mapping(1);
+  QGauss<dim>           quadrature(degree + 1);
+
+  DoFHandler<dim> dof_handler(tria);
+
+  unsigned int counter = 0;
+
+  for (auto &cell : dof_handler)
+    if (counter++ < tria.n_cells() / 2)
+      cell.set_active_fe_index(1);
+
+  dof_handler.distribute_dofs(fe);
+
+
+  typename MatrixFree<dim, double>::AdditionalData additional_data;
+  additional_data.mapping_update_flags_boundary_faces =
+    update_gradients | update_values;
+  additional_data.mapping_update_flags_inner_faces =
+    update_gradients | update_values;
+  additional_data.initialize_mapping = false;
+  additional_data.tasks_parallel_scheme =
+    MatrixFree<dim, double>::AdditionalData::none;
+
+  AffineConstraints<double> constraints;
+
+  MatrixFree<dim, double> matrix_free;
+  matrix_free.reinit(
+    mapping, dof_handler, constraints, quadrature, additional_data);
+
+  using VectorType = Vector<double>;
+
+  VectorType src, dst;
+  matrix_free.initialize_dof_vector(src);
+  matrix_free.initialize_dof_vector(dst);
+
+  std::vector<std::vector<CellId>>                          cells(2);
+  std::vector<std::vector<std::pair<CellId, unsigned int>>> boundary_faces(2);
+
+  std::vector<std::vector<std::vector<std::pair<CellId, CellId>>>> inner_faces(
+    2);
+  for (unsigned int i = 0; i < 2; ++i)
+    inner_faces[i].resize(2);
+
+  matrix_free.template loop<VectorType, VectorType>(
+    [&](const auto &data, auto &dst, const auto &src, const auto range) {
+      const auto i = data.get_cell_active_fe_index(range);
+
+      for (unsigned int cell = range.first; cell < range.second; ++cell)
+        for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell);
+             ++v)
+          cells[i].push_back(data.get_cell_iterator(cell, v)->id());
+    },
+    [&](const auto &data, auto &dst, const auto &src, const auto range) {
+      const auto i = data.get_face_active_fe_index(range, true);
+      const auto j = data.get_face_active_fe_index(range, false);
+
+      for (unsigned int face = range.first; face < range.second; ++face)
+        for (unsigned int v = 0; v < data.n_active_entries_per_face_batch(face);
+             ++v)
+          inner_faces[i][j].emplace_back(
+            data.get_face_iterator(face, v, true).first->id(),
+            data.get_face_iterator(face, v, false).first->id());
+    },
+    [&](const auto &data, auto &dst, const auto &src, const auto range) {
+      const auto i = data.get_face_active_fe_index(range);
+
+      for (unsigned int face = range.first; face < range.second; ++face)
+        for (unsigned int v = 0; v < data.n_active_entries_per_face_batch(face);
+             ++v)
+          boundary_faces[i].emplace_back(
+            data.get_face_iterator(face, v).first->id(),
+            data.get_face_iterator(face, v).second);
+    },
+    dst,
+    src);
+
+  deallog << "CELL categories:" << std::endl;
+  for (auto &i : cells)
+    {
+      std::sort(i.begin(), i.end());
+
+      for (const auto &j : i)
+        deallog << j << " ";
+      deallog << std::endl;
+    }
+
+  deallog << "BOUNDARY_FACE categories:" << std::endl;
+  for (auto &i : boundary_faces)
+    {
+      std::sort(i.begin(), i.end());
+      for (const auto &j : i)
+        deallog << j.first << "@" << j.second << "   ";
+      deallog << std::endl;
+    }
+
+  deallog << "INNER_FACE categories:" << std::endl;
+  for (auto &k : inner_faces)
+    {
+      for (auto &i : k)
+        {
+          std::sort(i.begin(), i.end());
+          for (const auto &j : i)
+            deallog << j.first << "@" << j.second << "   ";
+          deallog << std::endl;
+        }
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  test<2>();
+}
diff --git a/tests/simplex/matrix_free_range_iteration_01.output b/tests/simplex/matrix_free_range_iteration_01.output
new file mode 100644 (file)
index 0000000..3229c0c
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::CELL categories:
+DEAL::8_0: 9_0: 10_0: 11_0: 12_0: 13_0: 14_0: 15_0: 
+DEAL::0_0: 1_0: 2_0: 3_0: 4_0: 5_0: 6_0: 7_0: 
+DEAL::BOUNDARY_FACE categories:
+DEAL::8_0:@0   11_0:@1   12_0:@0   12_0:@3   13_0:@3   14_0:@3   15_0:@1   15_0:@3   
+DEAL::0_0:@0   0_0:@2   1_0:@2   2_0:@2   3_0:@1   3_0:@2   4_0:@0   7_0:@1   
+DEAL::INNER_FACE categories:
+DEAL::9_0:@8_0:   10_0:@9_0:   11_0:@10_0:   12_0:@8_0:   13_0:@9_0:   13_0:@12_0:   14_0:@10_0:   14_0:@13_0:   15_0:@11_0:   15_0:@14_0:   
+DEAL::
+DEAL::4_0:@8_0:   5_0:@9_0:   6_0:@10_0:   7_0:@11_0:   
+DEAL::1_0:@0_0:   2_0:@1_0:   3_0:@2_0:   4_0:@0_0:   5_0:@1_0:   5_0:@4_0:   6_0:@2_0:   6_0:@5_0:   7_0:@3_0:   7_0:@6_0:   

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.