add and test MatrixFree::get_cell_level_and_index()

diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h
index eec8213bea8208d943d5b633c8626aa52ba1a452..e70bdfbaaad3ecb140954e7fd37d37629587afa7 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.h
+++ b/include/deal.II/matrix_free/matrix_free.h
@@ -1455,6 +1455,15 @@ public:
                     const unsigned int vector_number,
                     const unsigned int dof_handler_index = 0) const;
 
+  /**
+   * This returns the level and index for the cell that would be
+   * returned by get_cell_iterator() for the same arguments @p
+   * macro_cell_number and @p vector_number.
+   */
+  std::pair<int, int>
+  get_cell_level_and_index(const unsigned int macro_cell_number,
+                           const unsigned int vector_number) const;
+
   /**
    * This returns the cell iterator in deal.II speak to a given cell in the
    * renumbering of this structure. This function returns an exception in case

diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h
index 23a95f10833c24969d6ba8f4d4ca2c2e1ceb1103..75ab6f44ab085d6b3b2d63f486e288f999e676b5 100644 (file)
--- a/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/include/deal.II/matrix_free/matrix_free.templates.h
@@ -188,6 +188,24 @@ MatrixFree<dim, Number, VectorizedArrayType>::get_cell_iterator(
 
 
 
+template <int dim, typename Number, typename VectorizedArrayType>
+std::pair<int, int>
+MatrixFree<dim, Number, VectorizedArrayType>::get_cell_level_and_index(
+  const unsigned int macro_cell_number,
+  const unsigned int vector_number) const
+{
+  const unsigned int vectorization_length =
+    VectorizedArrayType::n_array_elements;
+  AssertIndexRange(macro_cell_number, task_info.cell_partition_data.back());
+  AssertIndexRange(vector_number, n_components_filled(macro_cell_number));
+
+  std::pair<int, int> level_index_pair =
+    cell_level_index[macro_cell_number * vectorization_length + vector_number];
+  return level_index_pair;
+}
+
+
+
 template <int dim, typename Number, typename VectorizedArrayType>
 typename hp::DoFHandler<dim>::active_cell_iterator
 MatrixFree<dim, Number, VectorizedArrayType>::get_hp_cell_iterator(

diff --git a/tests/matrix_free/cell_level_and_index.cc b/tests/matrix_free/cell_level_and_index.cc
new file mode 100644 (file)
index 0000000..ac16715
--- /dev/null
+++ b/tests/matrix_free/cell_level_and_index.cc
@@ -0,0 +1,133 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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_cell_level_and_index()
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim, typename number>
+void
+compare_indices(const MatrixFree<dim, number> *mf_data)
+{
+  const unsigned int     n_batches  = mf_data->n_macro_cells();
+  constexpr unsigned int batch_size = VectorizedArray<number>::n_array_elements;
+  for (unsigned int batch_no = 0; batch_no < n_batches; ++batch_no)
+    {
+      const unsigned int n_lanes_filled =
+        mf_data->n_components_filled(batch_no);
+      for (unsigned int lane = 0; lane < n_lanes_filled; ++lane)
+        {
+          const auto cell = mf_data->get_cell_iterator(batch_no, lane);
+          const auto level_index_pair =
+            mf_data->get_cell_level_and_index(batch_no, lane);
+          std::cout << "(cell->level(), cell->index) = (" << cell->level()
+                    << ", " << cell->index() << ") "
+                    << "versus get_cell_level_and_index() = ("
+                    << level_index_pair.first << ", " << level_index_pair.second
+                    << ")\n";
+          AssertThrow((cell->level() == level_index_pair.first),
+                      ExcMessage("mismatching cell levels"));
+          AssertThrow((cell->index() == level_index_pair.second),
+                      ExcMessage("mismatching cell indices"));
+        }
+    }
+}
+
+
+
+template <int dim, int fe_degree, typename number = double>
+void
+test(const bool adaptive_ref = true)
+{
+  Triangulation<dim> tria(
+    Triangulation<dim>::limit_level_difference_at_vertices);
+  GridGenerator::hyper_cube(tria, 0, 1, true);
+  tria.refine_global(1);
+  if (adaptive_ref)
+    {
+      for (auto &cell : tria.active_cell_iterators())
+        if (cell->is_locally_owned())
+          if (cell->center().norm() < 0.5)
+            cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+      for (auto &cell : tria.active_cell_iterators())
+        if (cell->is_locally_owned())
+          if (cell->center()[0] < 0.2)
+            cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+  else
+    {
+      tria.refine_global(1);
+    }
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  dof.distribute_mg_dofs(fe);
+
+  AffineConstraints<double> constraints;
+  constraints.close();
+  const QGauss<1>                                  quad(fe_degree + 1);
+  typename MatrixFree<dim, number>::AdditionalData additional_data;
+  additional_data.tasks_parallel_scheme =
+    MatrixFree<dim, number>::AdditionalData::none;
+  auto mf_data = std::make_shared<MatrixFree<dim, number>>();
+
+  {
+    std::cout << "Compare active indices." << std::endl;
+    mf_data->reinit(dof, constraints, quad, additional_data);
+    compare_indices(mf_data.get());
+  }
+
+  {
+    std::cout << "Compare level indices." << std::endl;
+    const unsigned int level         = tria.n_global_levels() - 1;
+    additional_data.level_mg_handler = level;
+    mf_data->reinit(dof, constraints, quad, additional_data);
+    compare_indices(mf_data.get());
+  }
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  test<2, 1>();
+  test<2, 1, float>();
+  test<2, 1>(false);
+}

diff --git a/tests/matrix_free/cell_level_and_index.output b/tests/matrix_free/cell_level_and_index.output
new file mode 100644 (file)
index 0000000..089f418
--- /dev/null
+++ b/tests/matrix_free/cell_level_and_index.output
@@ -0,0 +1,86 @@
+Compare active indices.
+(cell->level(), cell->index) = (3, 0) versus get_cell_level_and_index() = (3, 0)
+(cell->level(), cell->index) = (3, 1) versus get_cell_level_and_index() = (3, 1)
+(cell->level(), cell->index) = (3, 2) versus get_cell_level_and_index() = (3, 2)
+(cell->level(), cell->index) = (3, 3) versus get_cell_level_and_index() = (3, 3)
+(cell->level(), cell->index) = (2, 1) versus get_cell_level_and_index() = (2, 1)
+(cell->level(), cell->index) = (3, 4) versus get_cell_level_and_index() = (3, 4)
+(cell->level(), cell->index) = (3, 5) versus get_cell_level_and_index() = (3, 5)
+(cell->level(), cell->index) = (3, 6) versus get_cell_level_and_index() = (3, 6)
+(cell->level(), cell->index) = (3, 7) versus get_cell_level_and_index() = (3, 7)
+(cell->level(), cell->index) = (2, 3) versus get_cell_level_and_index() = (2, 3)
+(cell->level(), cell->index) = (1, 1) versus get_cell_level_and_index() = (1, 1)
+(cell->level(), cell->index) = (2, 4) versus get_cell_level_and_index() = (2, 4)
+(cell->level(), cell->index) = (2, 5) versus get_cell_level_and_index() = (2, 5)
+(cell->level(), cell->index) = (2, 6) versus get_cell_level_and_index() = (2, 6)
+(cell->level(), cell->index) = (2, 7) versus get_cell_level_and_index() = (2, 7)
+(cell->level(), cell->index) = (1, 3) versus get_cell_level_and_index() = (1, 3)
+Compare level indices.
+(cell->level(), cell->index) = (3, 0) versus get_cell_level_and_index() = (3, 0)
+(cell->level(), cell->index) = (3, 1) versus get_cell_level_and_index() = (3, 1)
+(cell->level(), cell->index) = (3, 2) versus get_cell_level_and_index() = (3, 2)
+(cell->level(), cell->index) = (3, 3) versus get_cell_level_and_index() = (3, 3)
+(cell->level(), cell->index) = (3, 4) versus get_cell_level_and_index() = (3, 4)
+(cell->level(), cell->index) = (3, 5) versus get_cell_level_and_index() = (3, 5)
+(cell->level(), cell->index) = (3, 6) versus get_cell_level_and_index() = (3, 6)
+(cell->level(), cell->index) = (3, 7) versus get_cell_level_and_index() = (3, 7)
+Compare active indices.
+(cell->level(), cell->index) = (3, 0) versus get_cell_level_and_index() = (3, 0)
+(cell->level(), cell->index) = (3, 1) versus get_cell_level_and_index() = (3, 1)
+(cell->level(), cell->index) = (3, 2) versus get_cell_level_and_index() = (3, 2)
+(cell->level(), cell->index) = (3, 3) versus get_cell_level_and_index() = (3, 3)
+(cell->level(), cell->index) = (2, 1) versus get_cell_level_and_index() = (2, 1)
+(cell->level(), cell->index) = (3, 4) versus get_cell_level_and_index() = (3, 4)
+(cell->level(), cell->index) = (3, 5) versus get_cell_level_and_index() = (3, 5)
+(cell->level(), cell->index) = (3, 6) versus get_cell_level_and_index() = (3, 6)
+(cell->level(), cell->index) = (3, 7) versus get_cell_level_and_index() = (3, 7)
+(cell->level(), cell->index) = (2, 3) versus get_cell_level_and_index() = (2, 3)
+(cell->level(), cell->index) = (1, 1) versus get_cell_level_and_index() = (1, 1)
+(cell->level(), cell->index) = (2, 4) versus get_cell_level_and_index() = (2, 4)
+(cell->level(), cell->index) = (2, 5) versus get_cell_level_and_index() = (2, 5)
+(cell->level(), cell->index) = (2, 6) versus get_cell_level_and_index() = (2, 6)
+(cell->level(), cell->index) = (2, 7) versus get_cell_level_and_index() = (2, 7)
+(cell->level(), cell->index) = (1, 3) versus get_cell_level_and_index() = (1, 3)
+Compare level indices.
+(cell->level(), cell->index) = (3, 0) versus get_cell_level_and_index() = (3, 0)
+(cell->level(), cell->index) = (3, 1) versus get_cell_level_and_index() = (3, 1)
+(cell->level(), cell->index) = (3, 2) versus get_cell_level_and_index() = (3, 2)
+(cell->level(), cell->index) = (3, 3) versus get_cell_level_and_index() = (3, 3)
+(cell->level(), cell->index) = (3, 4) versus get_cell_level_and_index() = (3, 4)
+(cell->level(), cell->index) = (3, 5) versus get_cell_level_and_index() = (3, 5)
+(cell->level(), cell->index) = (3, 6) versus get_cell_level_and_index() = (3, 6)
+(cell->level(), cell->index) = (3, 7) versus get_cell_level_and_index() = (3, 7)
+Compare active indices.
+(cell->level(), cell->index) = (2, 0) versus get_cell_level_and_index() = (2, 0)
+(cell->level(), cell->index) = (2, 1) versus get_cell_level_and_index() = (2, 1)
+(cell->level(), cell->index) = (2, 2) versus get_cell_level_and_index() = (2, 2)
+(cell->level(), cell->index) = (2, 3) versus get_cell_level_and_index() = (2, 3)
+(cell->level(), cell->index) = (2, 4) versus get_cell_level_and_index() = (2, 4)
+(cell->level(), cell->index) = (2, 5) versus get_cell_level_and_index() = (2, 5)
+(cell->level(), cell->index) = (2, 6) versus get_cell_level_and_index() = (2, 6)
+(cell->level(), cell->index) = (2, 7) versus get_cell_level_and_index() = (2, 7)
+(cell->level(), cell->index) = (2, 8) versus get_cell_level_and_index() = (2, 8)
+(cell->level(), cell->index) = (2, 9) versus get_cell_level_and_index() = (2, 9)
+(cell->level(), cell->index) = (2, 10) versus get_cell_level_and_index() = (2, 10)
+(cell->level(), cell->index) = (2, 11) versus get_cell_level_and_index() = (2, 11)
+(cell->level(), cell->index) = (2, 12) versus get_cell_level_and_index() = (2, 12)
+(cell->level(), cell->index) = (2, 13) versus get_cell_level_and_index() = (2, 13)
+(cell->level(), cell->index) = (2, 14) versus get_cell_level_and_index() = (2, 14)
+(cell->level(), cell->index) = (2, 15) versus get_cell_level_and_index() = (2, 15)
+Compare level indices.
+(cell->level(), cell->index) = (2, 0) versus get_cell_level_and_index() = (2, 0)
+(cell->level(), cell->index) = (2, 1) versus get_cell_level_and_index() = (2, 1)
+(cell->level(), cell->index) = (2, 2) versus get_cell_level_and_index() = (2, 2)
+(cell->level(), cell->index) = (2, 3) versus get_cell_level_and_index() = (2, 3)
+(cell->level(), cell->index) = (2, 4) versus get_cell_level_and_index() = (2, 4)
+(cell->level(), cell->index) = (2, 5) versus get_cell_level_and_index() = (2, 5)
+(cell->level(), cell->index) = (2, 6) versus get_cell_level_and_index() = (2, 6)
+(cell->level(), cell->index) = (2, 7) versus get_cell_level_and_index() = (2, 7)
+(cell->level(), cell->index) = (2, 8) versus get_cell_level_and_index() = (2, 8)
+(cell->level(), cell->index) = (2, 9) versus get_cell_level_and_index() = (2, 9)
+(cell->level(), cell->index) = (2, 10) versus get_cell_level_and_index() = (2, 10)
+(cell->level(), cell->index) = (2, 11) versus get_cell_level_and_index() = (2, 11)
+(cell->level(), cell->index) = (2, 12) versus get_cell_level_and_index() = (2, 12)
+(cell->level(), cell->index) = (2, 13) versus get_cell_level_and_index() = (2, 13)
+(cell->level(), cell->index) = (2, 14) versus get_cell_level_and_index() = (2, 14)
+(cell->level(), cell->index) = (2, 15) versus get_cell_level_and_index() = (2, 15)