]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Access neighbors during ECL 9154/head
authorPeter Munch <peterrmuench@gmail.com>
Sun, 8 Dec 2019 13:59:54 +0000 (14:59 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Mon, 9 Dec 2019 15:34:25 +0000 (16:34 +0100)
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/matrix_free.h
include/deal.II/matrix_free/matrix_free.templates.h
tests/matrix_free/read_cell_data_01.cc [moved from tests/matrix_free/read_cell_data.cc with 100% similarity]
tests/matrix_free/read_cell_data_01.output [moved from tests/matrix_free/read_cell_data.output with 100% similarity]
tests/matrix_free/read_cell_data_02.cc [new file with mode: 0644]
tests/matrix_free/read_cell_data_02.with_mpi=true.with_p4est=true.mpirun=2.output [new file with mode: 0644]

index a654917d6752732d22ff78186898b1a9bca11c36..e8bb47f928d506f56d60a9c5ffcb7dfaa8d41d32 100644 (file)
@@ -611,6 +611,13 @@ public:
     const AlignedVector<std::array<T, VectorizedArrayType::n_array_elements>>
       &array) const;
 
+  /**
+   * Return the id of the cells this FEEvaluation or FEFaceEvaluation is
+   * associated with.
+   */
+  std::array<unsigned int, VectorizedArrayType::n_array_elements>
+  get_cell_ids() const;
+
   //@}
 
   /**
@@ -3555,6 +3562,82 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     return jacobian[q_index];
 }
 
+template <int dim,
+          int n_components_,
+          typename Number,
+          bool is_face,
+          typename VectorizedArrayType>
+std::array<unsigned int, VectorizedArrayType::n_array_elements>
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+  get_cell_ids() const
+{
+  const unsigned int v_len = VectorizedArrayType::n_array_elements;
+  std::array<unsigned int, VectorizedArrayType::n_array_elements> cells;
+
+  // initialize array
+  for (unsigned int i = 0; i < v_len; ++i)
+    cells[i] = numbers::invalid_unsigned_int;
+
+  if ((is_face == false) ||
+      (is_face &&
+       this->dof_access_index ==
+         internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+       this->is_interior_face))
+    {
+      // cell or interior face face (element-centric loop)
+      for (unsigned int i = 0; i < v_len; ++i)
+        cells[i] = cell * v_len + i;
+    }
+  else if (is_face &&
+           this->dof_access_index ==
+             internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+           this->is_interior_face == false)
+    {
+      // exterior face (element-centric loop): for this case, we need to
+      // look into the FaceInfo field that collects information from both
+      // sides of a face once for the global mesh, and pick the face id that
+      // is not the local one (cell_this).
+      for (unsigned int i = 0; i < v_len; i++)
+        {
+          // compute actual (non vectorized) cell ID
+          const unsigned int cell_this = this->cell * v_len + i;
+          // compute face ID
+          unsigned int face_index =
+            this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
+                                                                  this->face_no,
+                                                                  i);
+
+          if (face_index == numbers::invalid_unsigned_int)
+            continue; // invalid face ID: no neighbor on boundary
+
+          // get cell ID on both sides of face
+          auto cell_m = this->matrix_info->get_face_info(face_index / v_len)
+                          .cells_interior[face_index % v_len];
+          auto cell_p = this->matrix_info->get_face_info(face_index / v_len)
+                          .cells_exterior[face_index % v_len];
+
+          // compare the IDs with the given cell ID
+          if (cell_m == cell_this)
+            cells[i] = cell_p; // neighbor has the other ID
+          else if (cell_p == cell_this)
+            cells[i] = cell_m;
+        }
+    }
+  else if (is_face)
+    {
+      // face-centric faces
+      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)
+          cells[i] = cells_[i];
+    }
+
+  return cells;
+}
+
 
 
 template <int dim,
@@ -3569,21 +3652,17 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   Assert(matrix_info != nullptr, ExcNotImplemented());
   AssertDimension(array.size(),
                   matrix_info->get_task_info().cell_partition_data.back());
-  if (is_face)
-    {
-      VectorizedArrayType out = make_vectorized_array<Number>(Number(1.));
-      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];
+
+  // 1) collect ids of cell
+  const auto cells = this->get_cell_ids();
+
+  // 2) actually gather values
+  VectorizedArrayType out = make_vectorized_array<Number>(Number(1.));
+  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;
 }
 
 
@@ -3603,21 +3682,17 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   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];
+
+  // 1) collect ids of cell
+  const auto cells = this->get_cell_ids();
+
+  // 2) actually gather values
+  std::array<T, VectorizedArrayType::n_array_elements> out;
+  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;
 }
 
 
@@ -7313,10 +7388,6 @@ FEFaceEvaluation<dim,
   Assert(this->mapped_geometry == nullptr,
          ExcMessage("FEEvaluation was initialized without a matrix-free object."
                     " Integer indexing is not possible"));
-  Assert(this->is_interior_face == true,
-         ExcMessage(
-           "Cell-based FEFaceEvaluation::reinit only possible for the "
-           "interior face with second argument to constructor as true"));
   if (this->mapped_geometry != nullptr)
     return;
   Assert(this->matrix_info != nullptr, ExcNotInitialized());
index df97932f76de8a72248d885f2e38eae271ea20b9..4907ffe5494a3a84e5bed00212e2aa0b16e97238 100644 (file)
@@ -1778,6 +1778,15 @@ public:
     VectorizedArrayType::n_array_elements> &
   get_face_info(const unsigned int face_batch_number) const;
 
+
+  /**
+   * Return the table that translates a triple of the macro cell number,
+   * the index of a face within a cell and the index within the cell batch of
+   * vectorization into the index within the faces array.
+   */
+  const Table<3, unsigned int> &
+  get_cell_and_face_to_plain_faces() const;
+
   /**
    * Obtains a scratch data object for internal use. Make sure to release it
    * afterwards by passing the pointer you obtain from this object to the
@@ -2465,6 +2474,16 @@ MatrixFree<dim, Number, VectorizedArrayType>::get_face_info(
 
 
 
+template <int dim, typename Number, typename VectorizedArrayType>
+inline const Table<3, unsigned int> &
+MatrixFree<dim, Number, VectorizedArrayType>::get_cell_and_face_to_plain_faces()
+  const
+{
+  return face_info.cell_and_face_to_plain_faces;
+}
+
+
+
 template <int dim, typename Number, typename VectorizedArrayType>
 inline const Quadrature<dim> &
 MatrixFree<dim, Number, VectorizedArrayType>::get_quadrature(
index 8db758b14e9a8f2b8bf4bd1da4dddacb9f24323d..b52d71be148d10a21f691b13b5b34473a7906fdc 100644 (file)
@@ -1484,7 +1484,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
         true);
       face_info.cell_and_face_boundary_id.fill(numbers::invalid_boundary_id);
 
-      for (unsigned int f = 0; f < task_info.boundary_partition_data.back();
+      for (unsigned int f = 0; f < task_info.ghost_face_partition_data.back();
            ++f)
         for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements &&
                                  face_info.faces[f].cells_interior[v] !=
diff --git a/tests/matrix_free/read_cell_data_02.cc b/tests/matrix_free/read_cell_data_02.cc
new file mode 100644 (file)
index 0000000..d7aad75
--- /dev/null
@@ -0,0 +1,180 @@
+// ---------------------------------------------------------------------
+//
+// 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 each cell (incl. inner and outer faces,
+// which is needed for element-centric-loops)
+
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/tria.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)
+  {
+    parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+    GridGenerator::hyper_cube(tria);
+    tria.refine_global(2);
+
+    FE_DGQ<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> 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.mapping_update_flags_faces_by_cells = update_values;
+    additional_data.hold_all_faces_to_owned_cells       = true;
+
+    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.cell_loop(&Test::cell_wise_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_wise_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);
+    FEFaceEvaluation<dim, 1> fe_eval_m(data, true);
+    FEFaceEvaluation<dim, 1> fe_eval_p(data, false);
+
+    for (auto cell = pair.first; cell < pair.second; cell++)
+      {
+        for (auto lane = 0u; lane < data.n_active_entries_per_cell_batch(cell);
+             lane++)
+          {
+            fe_eval.reinit(cell);
+            const auto cell_data = fe_eval.read_cell_data(cell_ids);
+            Assert(cell_data[lane] == data.get_cell_iterator(cell, lane)->id(),
+                   ExcInternalError());
+            // deallog << "c " << cell_data[lane] << std::endl;
+
+            for (unsigned int face = 0;
+                 face < GeometryInfo<dim>::faces_per_cell;
+                 ++face)
+              {
+                fe_eval_m.reinit(cell, face);
+                fe_eval_p.reinit(cell, face);
+
+                const auto cell_data_m = fe_eval_m.read_cell_data(cell_ids);
+                Assert(cell_data_m[lane] ==
+                         data.get_cell_iterator(cell, lane)->id(),
+                       ExcInternalError());
+                // deallog << "m" << cell_data_m[lane] << std::endl;
+
+                const auto bids =
+                  data.get_faces_by_cells_boundary_id(cell, face);
+
+                if (bids[lane] != numbers::invalid_unsigned_int)
+                  continue;
+
+                const auto cell_data_p = fe_eval_p.read_cell_data(cell_ids);
+
+                Assert(
+                  cell_data_p[lane] ==
+                    data.get_cell_iterator(cell, lane)->neighbor(face)->id(),
+                  ExcInternalError());
+                // deallog << "p " << cell_data_p[lane] << std::endl;
+              }
+          }
+      }
+  }
+
+  AlignedVector<std::array<CellId, VectorizedArrayType::n_array_elements>>
+    cell_ids;
+};
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  MPILogInitAll all;
+
+  {
+    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_02.with_mpi=true.with_p4est=true.mpirun=2.output b/tests/matrix_free/read_cell_data_02.with_mpi=true.with_p4est=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..a924a83
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:0:2d::
+DEAL:0:3d::
+
+DEAL:1:2d::
+DEAL:1:3d::
+

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.