Introduce inverse of MF::get_cell_iterator()

diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h
index d8e519ab6d3d3e11788cbc6e97fe91c8ba891fcc..60334d1be256c94e2c798acaea45187f8af4d7f2 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.h
+++ b/include/deal.II/matrix_free/matrix_free.h
@@ -1656,6 +1656,16 @@ public:
   get_cell_level_and_index(const unsigned int cell_batch_index,
                            const unsigned int lane_index) const;
 
+  /**
+   * Get MatrixFree index associated to a deal.II @p cell. To get
+   * the actual cell batch index and lane, do the postprocessing
+   * `index / VectorizedArrayType::size()` and `index %
+   * VectorizedArrayType::size()`.
+   */
+  unsigned int
+  get_matrix_free_cell_index(
+    const typename Triangulation<dim>::cell_iterator &cell) const;
+
   /**
    * Return the cell iterator in deal.II speak to an interior/exterior cell of
    * a face in a pair of a face batch and lane index. The second element of
@@ -2043,6 +2053,11 @@ private:
    */
   std::vector<std::pair<unsigned int, unsigned int>> cell_level_index;
 
+  /**
+   * Conversion from deal.II index (active or level index) to MatrixFree index
+   * (inverse of cell_level_index).
+   */
+  std::vector<unsigned int> mf_cell_indices;
 
   /**
    * For discontinuous Galerkin, the cell_level_index includes cells that are

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 90417b5f0616753f9b03c851a760d89c8b558ca3..2e42dcbc39ce7f8e85edd2330a692782c5b22ee7 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -270,6 +270,19 @@ MatrixFree<dim, Number, VectorizedArrayType>::get_cell_level_and_index(
 
 
 
+template <int dim, typename Number, typename VectorizedArrayType>
+unsigned int
+MatrixFree<dim, Number, VectorizedArrayType>::get_matrix_free_cell_index(
+  const typename Triangulation<dim>::cell_iterator &cell) const
+{
+  return mf_cell_indices[(this->get_mg_level() ==
+                          numbers::invalid_unsigned_int) ?
+                           cell->active_cell_index() :
+                           cell->index()];
+}
+
+
+
 template <int dim, typename Number, typename VectorizedArrayType>
 std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>
 MatrixFree<dim, Number, VectorizedArrayType>::get_face_iterator(
@@ -802,6 +815,29 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
 
       mapping_is_initialized = true;
     }
+
+  // set up map: deal.II index -> MatrixFree index
+  {
+    const auto &tria     = dof_handler[0]->get_triangulation();
+    const auto  mg_level = this->get_mg_level();
+
+    mf_cell_indices.resize((mg_level == numbers::invalid_unsigned_int) ?
+                             tria.n_active_cells() :
+                             tria.n_cells(mg_level),
+                           numbers::invalid_unsigned_int);
+
+    for (unsigned int cell = 0;
+         cell < n_cell_batches() + n_ghost_cell_batches();
+         ++cell)
+      for (unsigned int v = 0; v < n_active_entries_per_cell_batch(cell); ++v)
+        {
+          const auto tria_cell = get_cell_iterator(cell, v);
+          mf_cell_indices[(mg_level == numbers::invalid_unsigned_int) ?
+                            tria_cell->active_cell_index() :
+                            tria_cell->index()] =
+            cell * VectorizedArrayType::size() + v;
+        }
+  }
 }
 
 

diff --git a/tests/matrix_free/iterators.cc b/tests/matrix_free/iterators_01.cc
similarity index 100%
rename from tests/matrix_free/iterators.cc
rename to tests/matrix_free/iterators_01.cc

diff --git a/tests/matrix_free/iterators.output b/tests/matrix_free/iterators_01.output
similarity index 100%
rename from tests/matrix_free/iterators.output
rename to tests/matrix_free/iterators_01.output

diff --git a/tests/matrix_free/iterators_02.cc b/tests/matrix_free/iterators_02.cc
new file mode 100644 (file)
index 0000000..ead6066
--- /dev/null
+++ b/tests/matrix_free/iterators_02.cc
@@ -0,0 +1,96 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Test MatrixFree::get_matrix_free_cell_index() and
+// MatrixFree::get_cell_iterator(), whether they indeed return
+// inverse information.
+
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <set>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test()
+{
+  const unsigned int fe_degree = 1;
+
+  using Number              = double;
+  using VectorizedArrayType = VectorizedArray<Number>;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(3);
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+
+  QGauss<1> quad(fe_degree + 1);
+
+  AffineConstraints<Number> constraints;
+
+  MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+  matrix_free.reinit(mapping, dof_handler, constraints, quad);
+
+  for (const auto &cell : tria.active_cell_iterators())
+    {
+      const auto mf_index = matrix_free.get_matrix_free_cell_index(cell);
+
+      AssertThrow(cell == matrix_free.get_cell_iterator(
+                            mf_index / VectorizedArrayType::size(),
+                            mf_index % VectorizedArrayType::size()),
+                  ExcInternalError());
+    }
+
+  deallog << "OK!" << std::endl;
+
+  for (unsigned int i = 0; i < matrix_free.n_cell_batches(); ++i)
+    for (unsigned int v = 0; v < matrix_free.n_active_entries_per_cell_batch(i);
+         ++v)
+      AssertDimension(matrix_free.get_matrix_free_cell_index(
+                        matrix_free.get_cell_iterator(i, v)),
+                      i * VectorizedArrayType::size() + v);
+
+  deallog << "OK!" << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+
+  test<2>();
+}

diff --git a/tests/matrix_free/iterators_02.output b/tests/matrix_free/iterators_02.output
new file mode 100644 (file)
index 0000000..043986a
--- /dev/null
+++ b/tests/matrix_free/iterators_02.output
@@ -0,0 +1,3 @@
+
+DEAL::OK!
+DEAL::OK!