From: Peter Munch <peterrmuench@gmail.com>
Date: Sun, 23 Jan 2022 17:50:59 +0000 (+0100)
Subject: Introduce FEEval::*_face_data()
X-Git-Tag: v9.4.0-rc1~546^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e5148c2ee603e7556d1935efc104fbc5cfcfb928;p=dealii.git

Introduce FEEval::*_face_data()
---

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 1e2c7760a2..030932d6c6 100644
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -579,44 +579,6 @@ public:
 
   //@}
 
-  /**
-   * Provides a unified interface to access data in a vector of
-   * VectorizedArray fields of length MatrixFree::n_cell_batches() +
-   * MatrixFree::n_ghost_cell_batches() for both cells (plain read) and faces
-   * (indirect addressing).
-   */
-  VectorizedArrayType
-  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const;
-
-  /**
-   * Provides a unified interface to set data in a vector of
-   * VectorizedArray fields of length MatrixFree::n_cell_batches() +
-   * MatrixFree::n_ghost_cell_batches() for both cells (plain read) and faces
-   * (indirect addressing).
-   */
-  void
-  set_cell_data(AlignedVector<VectorizedArrayType> &array,
-                const VectorizedArrayType &         value) const;
-
-  /**
-   * The same as above, just for std::array of length of VectorizedArrayType for
-   * arbitrary data type.
-   */
-  template <typename T>
-  std::array<T, VectorizedArrayType::size()>
-  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
-                   &array) const;
-
-  /**
-   * The same as above, just for std::array of length of VectorizedArrayType for
-   * arbitrary data type.
-   */
-  template <typename T>
-  void
-  set_cell_data(
-    AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
-    const std::array<T, VectorizedArrayType::size()> &         value) const;
-
   /**
    * Return the underlying MatrixFree object.
    */
@@ -3144,124 +3106,6 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 }
 
 
-namespace internal
-{
-  template <int dim,
-            int n_components,
-            typename Number,
-            bool is_face,
-            typename VectorizedArrayType,
-            typename VectorizedArrayType2,
-            typename GlobalVectorType,
-            typename FU>
-  inline void
-  process_cell_data(
-    const FEEvaluationBase<dim,
-                           n_components,
-                           Number,
-                           is_face,
-                           VectorizedArrayType> &       fe_eval,
-    const MatrixFree<dim, Number, VectorizedArrayType> *matrix_free,
-    GlobalVectorType &                                  array,
-    VectorizedArrayType2 &                              out,
-    const FU &                                          fu)
-  {
-    (void)matrix_free;
-    Assert(matrix_free != nullptr, ExcNotImplemented());
-    AssertDimension(array.size(),
-                    matrix_free->get_task_info().cell_partition_data.back());
-
-    // 1) collect ids of cell
-    const auto cells = fe_eval.get_cell_ids();
-
-    // 2) actually gather values
-    for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
-      if (cells[i] != numbers::invalid_unsigned_int)
-        fu(out[i],
-           array[cells[i] / VectorizedArrayType::size()]
-                [cells[i] % VectorizedArrayType::size()]);
-  }
-} // namespace internal
-
-
-
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-inline VectorizedArrayType
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const
-{
-  VectorizedArrayType out = Number(1.);
-  internal::process_cell_data(
-    *this, this->matrix_free, array, out, [](auto &local, const auto &global) {
-      local = global;
-    });
-  return out;
-}
-
-
-
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-inline void
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  set_cell_data(AlignedVector<VectorizedArrayType> &array,
-                const VectorizedArrayType &         in) const
-{
-  internal::process_cell_data(
-    *this, this->matrix_free, array, in, [](const auto &local, auto &global) {
-      global = local;
-    });
-}
-
-
-
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-template <typename T>
-inline std::array<T, VectorizedArrayType::size()>
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
-                   &array) const
-{
-  std::array<T, VectorizedArrayType::size()> out;
-  internal::process_cell_data(
-    *this, this->matrix_free, array, out, [](auto &local, const auto &global) {
-      local = global;
-    });
-  return out;
-}
-
-
-
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-template <typename T>
-inline void
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  set_cell_data(
-    AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
-    const std::array<T, VectorizedArrayType::size()> &         in) const
-{
-  internal::process_cell_data(
-    *this, this->matrix_free, array, in, [](const auto &local, auto &global) {
-      global = local;
-    });
-}
-
-
 
 namespace internal
 {
@@ -7583,9 +7427,9 @@ FEFaceEvaluation<dim,
   unsigned int i = 0;
   for (; i < this->matrix_free->n_active_entries_per_face_batch(this->cell);
        ++i)
-    this->cell_or_face_ids[i] = face_index * VectorizedArrayType::size() + i;
+    this->face_ids[i] = face_index * VectorizedArrayType::size() + i;
   for (; i < VectorizedArrayType::size(); ++i)
-    this->cell_or_face_ids[i] = numbers::invalid_unsigned_int;
+    this->face_ids[i] = numbers::invalid_unsigned_int;
 
   this->cell_type = this->matrix_free->get_mapping_info().face_type[face_index];
   const unsigned int offsets =
@@ -7674,6 +7518,8 @@ FEFaceEvaluation<dim,
                                                                   face_number,
                                                                   i);
 
+          this->face_ids[i] = face_index;
+
           if (face_index == numbers::invalid_unsigned_int)
             {
               this->cell_ids[i]          = numbers::invalid_unsigned_int;
@@ -7722,6 +7568,11 @@ FEFaceEvaluation<dim,
       this->face_numbers[0]      = face_number;
       for (unsigned int i = 0; i < n_lanes; ++i)
         this->cell_ids[i] = cell_index * n_lanes + i;
+      for (unsigned int i = 0; i < n_lanes; ++i)
+        this->face_ids[i] =
+          this->matrix_free->get_cell_and_face_to_plain_faces()(cell_index,
+                                                                face_number,
+                                                                i);
     }
 
   const unsigned int offsets =
diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h
index 7de101a44f..2661c074d5 100644
--- a/include/deal.II/matrix_free/fe_evaluation_data.h
+++ b/include/deal.II/matrix_free/fe_evaluation_data.h
@@ -485,6 +485,20 @@ public:
     return cell_ids;
   }
 
+  /**
+   * Return the id of the faces this FEFaceEvaluation is
+   * associated with.
+   */
+  const std::array<unsigned int, n_lanes> &
+  get_face_ids() const
+  {
+// implemented inline to avoid compilation problems on Windows
+#ifdef DEBUG
+    Assert(is_reinitialized && is_face, ExcNotInitialized());
+#endif
+    return face_ids;
+  }
+
   /**
    * Return the id of the cell/face batch this FEEvaluation/FEFaceEvaluation is
    * associated with.
@@ -516,11 +530,94 @@ public:
                       internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
       return cell_ids;
     else
-      return cell_or_face_ids;
+      return face_ids;
   }
 
   //@}
 
+  /**
+   * @name 3: Functions to access cell- and face-data vectors.
+   */
+  //@{
+
+  /**
+   * Provides a unified interface to access data in a vector of
+   * VectorizedArray fields of length MatrixFree::n_cell_batches() +
+   * MatrixFree::n_ghost_cell_batches() for cell data. It is implemented
+   * both for cells and faces (access data to the associated cell).
+   */
+  Number
+  read_cell_data(const AlignedVector<Number> &array) const;
+
+  /**
+   * Provides a unified interface to set data in a vector of
+   * VectorizedArray fields of length MatrixFree::n_cell_batches() +
+   * MatrixFree::n_ghost_cell_batches() for cell data. It is implemented
+   * both for cells and faces (access data to the associated cell).
+   */
+  void
+  set_cell_data(AlignedVector<Number> &array, const Number &value) const;
+
+  /**
+   * The same as above, just for std::array of length of Number for
+   * arbitrary data type.
+   */
+  template <typename T>
+  std::array<T, Number::size()>
+  read_cell_data(
+    const AlignedVector<std::array<T, Number::size()>> &array) const;
+
+  /**
+   * The same as above, just for std::array of length of Number for
+   * arbitrary data type.
+   */
+  template <typename T>
+  void
+  set_cell_data(AlignedVector<std::array<T, Number::size()>> &array,
+                const std::array<T, Number::size()> &         value) const;
+
+  /**
+   * Provides a unified interface to access data in a vector of
+   * VectorizedArray fields of length MatrixFree::n_inner_face_batches() +
+   * MatrixFree::n_boundary_face_batches() +
+   * MatrixFree::n_ghost_inner_face_batches() for face data.
+   *
+   * @note Only implemented for faces.
+   */
+  Number
+  read_face_data(const AlignedVector<Number> &array) const;
+
+  /**
+   * Provides a unified interface to set data in a vector of
+   * VectorizedArray fields of length MatrixFree::n_inner_face_batches() +
+   * MatrixFree::n_boundary_face_batches() +
+   * MatrixFree::n_ghost_inner_face_batches() for face data.
+   *
+   * @note Only implemented for faces.
+   */
+  void
+  set_face_data(AlignedVector<Number> &array, const Number &value) const;
+
+  /**
+   * The same as above, just for std::array of length of Number for
+   * arbitrary data type.
+   */
+  template <typename T>
+  std::array<T, Number::size()>
+  read_face_data(
+    const AlignedVector<std::array<T, Number::size()>> &array) const;
+
+  /**
+   * The same as above, just for std::array of length of Number for
+   * arbitrary data type.
+   */
+  template <typename T>
+  void
+  set_face_data(AlignedVector<std::array<T, Number::size()>> &array,
+                const std::array<T, Number::size()> &         value) const;
+
+  //@}
+
   /**
    * This data structure is used for the initialization by the derived
    * FEEvaluation classes.
@@ -855,7 +952,7 @@ protected:
    * Stores the (non-vectorized) id of the cells or faces this object is
    * initialized to. Relevant for ECL.
    */
-  std::array<unsigned int, n_lanes> cell_or_face_ids;
+  std::array<unsigned int, n_lanes> face_ids;
 
   /**
    * Geometry data that can be generated FEValues on the fly with the
@@ -1485,6 +1582,165 @@ FEEvaluationData<dim, Number, is_face>::get_is_interior_face() const
 
 
 
+namespace internal
+{
+  template <std::size_t N,
+            typename VectorizedArrayType2,
+            typename GlobalVectorType,
+            typename FU>
+  inline void
+  process_cell_or_face_data(const std::array<unsigned int, N> indices,
+                            GlobalVectorType &                array,
+                            VectorizedArrayType2 &            out,
+                            const FU &                        fu)
+  {
+    for (unsigned int i = 0; i < N; ++i)
+      if (indices[i] != numbers::invalid_unsigned_int)
+        {
+          AssertIndexRange(indices[i] / N, array.size());
+          fu(out[i], array[indices[i] / N][indices[i] % N]);
+        }
+  }
+} // namespace internal
+
+
+
+template <int dim, typename Number, bool is_face>
+inline Number
+FEEvaluationData<dim, Number, is_face>::read_cell_data(
+  const AlignedVector<Number> &array) const
+{
+  Number out = Number(1.);
+  internal::process_cell_or_face_data(this->get_cell_ids(),
+                                      array,
+                                      out,
+                                      [](auto &local, const auto &global) {
+                                        local = global;
+                                      });
+  return out;
+}
+
+
+
+template <int dim, typename Number, bool is_face>
+inline void
+FEEvaluationData<dim, Number, is_face>::set_cell_data(
+  AlignedVector<Number> &array,
+  const Number &         in) const
+{
+  internal::process_cell_or_face_data(this->get_cell_ids(),
+                                      array,
+                                      in,
+                                      [](const auto &local, auto &global) {
+                                        global = local;
+                                      });
+}
+
+
+
+template <int dim, typename Number, bool is_face>
+template <typename T>
+inline std::array<T, Number::size()>
+FEEvaluationData<dim, Number, is_face>::read_cell_data(
+  const AlignedVector<std::array<T, Number::size()>> &array) const
+{
+  std::array<T, Number::size()> out;
+  internal::process_cell_or_face_data(this->get_cell_ids(),
+                                      array,
+                                      out,
+                                      [](auto &local, const auto &global) {
+                                        local = global;
+                                      });
+  return out;
+}
+
+
+
+template <int dim, typename Number, bool is_face>
+template <typename T>
+inline void
+FEEvaluationData<dim, Number, is_face>::set_cell_data(
+  AlignedVector<std::array<T, Number::size()>> &array,
+  const std::array<T, Number::size()> &         in) const
+{
+  internal::process_cell_or_face_data(this->get_cell_ids(),
+                                      array,
+                                      in,
+                                      [](const auto &local, auto &global) {
+                                        global = local;
+                                      });
+}
+
+
+
+template <int dim, typename Number, bool is_face>
+inline Number
+FEEvaluationData<dim, Number, is_face>::read_face_data(
+  const AlignedVector<Number> &array) const
+{
+  Number out = Number(1.);
+  internal::process_cell_or_face_data(this->get_face_ids(),
+                                      array,
+                                      out,
+                                      [](auto &local, const auto &global) {
+                                        local = global;
+                                      });
+  return out;
+}
+
+
+
+template <int dim, typename Number, bool is_face>
+inline void
+FEEvaluationData<dim, Number, is_face>::set_face_data(
+  AlignedVector<Number> &array,
+  const Number &         in) const
+{
+  internal::process_cell_or_face_data(this->get_face_ids(),
+                                      array,
+                                      in,
+                                      [](const auto &local, auto &global) {
+                                        global = local;
+                                      });
+}
+
+
+
+template <int dim, typename Number, bool is_face>
+template <typename T>
+inline std::array<T, Number::size()>
+FEEvaluationData<dim, Number, is_face>::read_face_data(
+  const AlignedVector<std::array<T, Number::size()>> &array) const
+{
+  std::array<T, Number::size()> out;
+  internal::process_cell_or_face_data(this->get_face_ids(),
+                                      array,
+                                      out,
+                                      [](auto &local, const auto &global) {
+                                        local = global;
+                                      });
+  return out;
+}
+
+
+
+template <int dim, typename Number, bool is_face>
+template <typename T>
+inline void
+FEEvaluationData<dim, Number, is_face>::set_face_data(
+  AlignedVector<std::array<T, Number::size()>> &array,
+  const std::array<T, Number::size()> &         in) const
+{
+  internal::process_cell_or_face_data(this->get_face_ids(),
+                                      array,
+                                      in,
+                                      [](const auto &local, auto &global) {
+                                        global = local;
+                                      });
+}
+
+
+
 #endif // ifndef DOXYGEN