Add new read_cell_data to FEEvaluationBase for arbitrary data types

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 9716922d9d11c8a5e0268633fca66072f6a045cf..0bb4baa57f3b2139757ab39f4b42ad1fe93a589a 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -601,6 +601,16 @@ public:
   VectorizedArrayType
   read_cell_data(const AlignedVector<VectorizedArrayType> &array) const;
 
+  /**
+   * The same as above, just for std::array of length of VectorizedArrayType for
+   * arbitrary data type.
+   */
+  template <typename T>
+  std::array<T, VectorizedArrayType::n_array_elements>
+  read_cell_data(
+    const AlignedVector<std::array<T, VectorizedArrayType::n_array_elements>>
+      &array) const;
+
   //@}
 
   /**
@@ -3578,6 +3588,40 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
 
 
+template <int dim,
+          int n_components_,
+          typename Number,
+          bool is_face,
+          typename VectorizedArrayType>
+template <typename T>
+inline std::array<T, VectorizedArrayType::n_array_elements>
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+  read_cell_data(
+    const AlignedVector<std::array<T, VectorizedArrayType::n_array_elements>>
+      &array) const
+{
+  Assert(matrix_info != nullptr, ExcNotImplemented());
+  AssertDimension(array.size(),
+                  matrix_info->get_task_info().cell_partition_data.back());
+  if (is_face)
+    {
+      std::array<T, VectorizedArrayType::n_array_elements> out;
+      const unsigned int *                                 cells =
+        is_interior_face ?
+          &this->matrix_info->get_face_info(cell).cells_interior[0] :
+          &this->matrix_info->get_face_info(cell).cells_exterior[0];
+      for (unsigned int i = 0; i < VectorizedArrayType::n_array_elements; ++i)
+        if (cells[i] != numbers::invalid_unsigned_int)
+          out[i] = array[cells[i] / VectorizedArrayType::n_array_elements]
+                        [cells[i] % VectorizedArrayType::n_array_elements];
+      return out;
+    }
+  else
+    return array[cell];
+}
+
+
+
 namespace internal
 {
   // allows to select between block vectors and non-block vectors, which

diff --git a/tests/matrix_free/read_cell_data.cc b/tests/matrix_free/read_cell_data.cc
new file mode 100644 (file)
index 0000000..4152562
--- /dev/null
+++ b/tests/matrix_free/read_cell_data.cc
@@ -0,0 +1,204 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// tests matrix-free read_cell_data for locally refined mesh
+
+
+#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_dgq.h>
+#include <deal.II/fe/mapping_q.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 <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim,
+          typename Number              = double,
+          typename VectorizedArrayType = VectorizedArray<Number>>
+class Test
+{
+public:
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  void
+  run(unsigned int fe_degree)
+  {
+    Triangulation<dim> tria(
+      Triangulation<dim>::limit_level_difference_at_vertices);
+    GridGenerator::hyper_cube(tria);
+
+    for (unsigned int i = 0; i < 3; i++)
+      {
+        tria.begin_active(i)->set_refine_flag();
+        tria.execute_coarsening_and_refinement();
+      }
+
+    FE_DGQ<dim>     fe(fe_degree);
+    DoFHandler<dim> dof_handler(tria);
+    dof_handler.distribute_dofs(fe);
+    dof_handler.distribute_mg_dofs();
+
+    MappingQ<dim> mapping(1);
+
+    QGauss<1> quad(fe_degree + 1);
+
+
+
+    unsigned int level = numbers::invalid_unsigned_int;
+    {
+      AffineConstraints<Number> constraint;
+
+      typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData
+        additional_data;
+      additional_data.mapping_update_flags                = update_values;
+      additional_data.mapping_update_flags_inner_faces    = update_values;
+      additional_data.mapping_update_flags_boundary_faces = update_values;
+      additional_data.mg_level                            = level;
+
+      MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+      matrix_free.reinit(
+        mapping, dof_handler, constraint, quad, additional_data);
+
+      VectorType src, dst;
+
+      matrix_free.initialize_dof_vector(src);
+      matrix_free.initialize_dof_vector(dst);
+
+      this->setup_vector(matrix_free);
+
+      matrix_free.loop(&Test::cell_operation,
+                       &Test::face_operation,
+                       &Test::boundary_operation,
+                       this,
+                       dst,
+                       src);
+
+      deallog << std::endl;
+    }
+  }
+
+private:
+  void
+  setup_vector(const MatrixFree<dim, Number, VectorizedArrayType> &data)
+  {
+    unsigned int n_cells = data.n_cell_batches() + data.n_ghost_cell_batches();
+    cell_ids.resize(n_cells);
+
+    for (unsigned int cell = 0; cell < n_cells; cell++)
+      for (auto lane = 0u; lane < data.n_active_entries_per_cell_batch(cell);
+           lane++)
+        cell_ids[cell][lane] = data.get_cell_iterator(cell, lane)->id();
+  }
+
+  void
+  cell_operation(const MatrixFree<dim, Number, VectorizedArrayType> &data,
+                 VectorType &,
+                 const VectorType &,
+                 const std::pair<unsigned int, unsigned int> &pair) const
+  {
+    FEEvaluation<dim, 1> fe_eval(data);
+    for (auto cell = pair.first; cell < pair.second; cell++)
+      {
+        fe_eval.reinit(cell);
+        const auto cell_data = fe_eval.read_cell_data(cell_ids);
+        for (auto lane = 0u; lane < data.n_active_entries_per_cell_batch(cell);
+             lane++)
+          Assert(cell_data[lane] == data.get_cell_iterator(cell, lane)->id(),
+                 ExcInternalError());
+      }
+  }
+
+  void
+  face_operation(const MatrixFree<dim, Number, VectorizedArrayType> &data,
+                 VectorType &,
+                 const VectorType &,
+                 const std::pair<unsigned int, unsigned int> &pair) const
+  {
+    FEFaceEvaluation<dim, 1> fe_eval_m(data, true);
+    FEFaceEvaluation<dim, 1> fe_eval_p(data, false);
+    for (auto face = pair.first; face < pair.second; face++)
+      {
+        fe_eval_m.reinit(face);
+        fe_eval_p.reinit(face);
+        const auto cell_data_m = fe_eval_m.read_cell_data(cell_ids);
+        const auto cell_data_p = fe_eval_p.read_cell_data(cell_ids);
+        for (auto lane = 0u; lane < data.n_active_entries_per_face_batch(face);
+             lane++)
+          {
+            Assert(cell_data_m[lane] ==
+                     data.get_face_iterator(face, lane, true).first->id(),
+                   ExcInternalError());
+            Assert(cell_data_p[lane] ==
+                     data.get_face_iterator(face, lane, false).first->id(),
+                   ExcInternalError());
+          }
+      }
+  }
+
+
+  void
+  boundary_operation(const MatrixFree<dim, Number, VectorizedArrayType> &data,
+                     VectorType &,
+                     const VectorType &,
+                     const std::pair<unsigned int, unsigned int> &pair) const
+  {
+    FEFaceEvaluation<dim, 1> fe_eval(data, true);
+    for (auto face = pair.first; face < pair.second; face++)
+      {
+        fe_eval.reinit(face);
+        const auto cell_data = fe_eval.read_cell_data(cell_ids);
+        for (auto lane = 0u; lane < data.n_active_entries_per_face_batch(face);
+             lane++)
+          Assert(cell_data[lane] ==
+                   data.get_face_iterator(face, lane, true).first->id(),
+                 ExcInternalError());
+      }
+  }
+
+  AlignedVector<std::array<CellId, VectorizedArrayType::n_array_elements>>
+    cell_ids;
+};
+
+
+int
+main()
+{
+  initlog();
+  {
+    deallog.push("2d");
+    Test<2> runner;
+    runner.run(1);
+    deallog.pop();
+  }
+  {
+    deallog.push("3d");
+    Test<3> runner;
+    runner.run(1);
+    deallog.pop();
+  }
+}

diff --git a/tests/matrix_free/read_cell_data.output b/tests/matrix_free/read_cell_data.output
new file mode 100644 (file)
index 0000000..6ec9f05
--- /dev/null
+++ b/tests/matrix_free/read_cell_data.output
@@ -0,0 +1,3 @@
+
+DEAL:2d::
+DEAL:3d::