From 23893199eb9ca63c1681fad7d30b9707c5157067 Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Fri, 27 Nov 2020 15:41:40 +0100
Subject: [PATCH] Add new constructors to FE(Face)Evaluation

---
 include/deal.II/matrix_free/fe_evaluation.h | 112 ++++++++++++++++++--
 tests/simplex/matrix_free_02.cc             |   6 +-
 2 files changed, 104 insertions(+), 14 deletions(-)

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 8583254d03..2018a77a27 100644
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -2556,6 +2556,19 @@ public:
     const unsigned int active_fe_index   = numbers::invalid_unsigned_int,
     const unsigned int active_quad_index = numbers::invalid_unsigned_int);
 
+  /**
+   * Constructor. Takes all data stored in MatrixFree for a given cell range,
+   * which allows to automatically identify the active_fe_index and
+   * active_quad_index in case of a p-adaptive strategy.
+   *
+   * The rest of the arguments are the same as in the constructor above.
+   */
+  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+               const std::pair<unsigned int, unsigned int> &       range,
+               const unsigned int                                  dof_no  = 0,
+               const unsigned int                                  quad_no = 0,
+               const unsigned int first_selected_component                 = 0);
+
   /**
    * Constructor that comes with reduced functionality and works similar as
    * FEValues. The arguments are similar to the ones passed to the constructor
@@ -3040,6 +3053,21 @@ public:
     const unsigned int active_fe_index   = numbers::invalid_unsigned_int,
     const unsigned int active_quad_index = numbers::invalid_unsigned_int);
 
+  /**
+   * Constructor. Takes all data stored in MatrixFree for a given face range,
+   * which allows to automatically identify the active_fe_index and
+   * active_quad_index in case of a p-adaptive strategy.
+   *
+   * The rest of the arguments are the same as in the constructor above.
+   */
+  FEFaceEvaluation(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const std::pair<unsigned int, unsigned int> &       range,
+    const bool                                          is_interior_face = true,
+    const unsigned int                                  dof_no           = 0,
+    const unsigned int                                  quad_no          = 0,
+    const unsigned int first_selected_component                          = 0);
+
   /**
    * Initializes the operation pointer to the current face. This method is the
    * default choice for face integration as the data stored in MappingInfo is
@@ -3295,16 +3323,24 @@ inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
                       (active_fe_index_in != numbers::invalid_unsigned_int ?
                          active_fe_index_in :
                          0))
-  , active_quad_index(fe_degree != numbers::invalid_unsigned_int ?
-                        (is_face ? data_in.get_mapping_info()
-                                     .face_data[quad_no_in]
-                                     .quad_index_from_n_q_points(n_q_points) :
-                                   data_in.get_mapping_info()
-                                     .cell_data[quad_no_in]
-                                     .quad_index_from_n_q_points(n_q_points)) :
-                        (active_quad_index_in != numbers::invalid_unsigned_int ?
-                           active_quad_index_in :
-                           0))
+  , active_quad_index(
+      fe_degree != numbers::invalid_unsigned_int ?
+        (is_face ? data_in.get_mapping_info()
+                     .face_data[quad_no_in]
+                     .quad_index_from_n_q_points(n_q_points) :
+                   data_in.get_mapping_info()
+                     .cell_data[quad_no_in]
+                     .quad_index_from_n_q_points(n_q_points)) :
+        (active_quad_index_in != numbers::invalid_unsigned_int ?
+           active_quad_index_in :
+           std::min<unsigned int>(active_fe_index,
+                                  (is_face ? data_in.get_mapping_info()
+                                               .face_data[quad_no_in]
+                                               .descriptor.size() :
+                                             data_in.get_mapping_info()
+                                               .cell_data[quad_no_in]
+                                               .descriptor.size()) -
+                                    1)))
   , n_quadrature_points(
       fe_degree != numbers::invalid_unsigned_int ?
         n_q_points :
@@ -7210,6 +7246,32 @@ inline FEEvaluation<dim,
 
 
 
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d,
+          int n_components_,
+          typename Number,
+          typename VectorizedArrayType>
+inline FEEvaluation<dim,
+                    fe_degree,
+                    n_q_points_1d,
+                    n_components_,
+                    Number,
+                    VectorizedArrayType>::
+  FEEvaluation(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+               const std::pair<unsigned int, unsigned int> &       range,
+               const unsigned int                                  dof_no,
+               const unsigned int                                  quad_no,
+               const unsigned int first_selected_component)
+  : FEEvaluation(matrix_free,
+                 dof_no,
+                 quad_no,
+                 first_selected_component,
+                 matrix_free.get_cell_active_fe_index(range))
+{}
+
+
+
 template <int dim,
           int fe_degree,
           int n_q_points_1d,
@@ -8321,6 +8383,36 @@ inline FEFaceEvaluation<dim,
 
 
 
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d,
+          int n_components_,
+          typename Number,
+          typename VectorizedArrayType>
+inline FEFaceEvaluation<dim,
+                        fe_degree,
+                        n_q_points_1d,
+                        n_components_,
+                        Number,
+                        VectorizedArrayType>::
+  FEFaceEvaluation(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const std::pair<unsigned int, unsigned int> &       range,
+    const bool                                          is_interior_face,
+    const unsigned int                                  dof_no,
+    const unsigned int                                  quad_no,
+    const unsigned int first_selected_component)
+  : FEFaceEvaluation(matrix_free,
+                     is_interior_face,
+                     dof_no,
+                     quad_no,
+                     first_selected_component,
+                     matrix_free.get_face_active_fe_index(range,
+                                                          is_interior_face))
+{}
+
+
+
 template <int dim,
           int fe_degree,
           int n_q_points_1d,
diff --git a/tests/simplex/matrix_free_02.cc b/tests/simplex/matrix_free_02.cc
index d884c35976..b510772772 100644
--- a/tests/simplex/matrix_free_02.cc
+++ b/tests/simplex/matrix_free_02.cc
@@ -183,8 +183,7 @@ public:
 
     matrix_free.template cell_loop<VectorType, int>(
       [&](const auto &data, auto &dst, const auto &, const auto range) {
-        const auto       i = data.get_cell_active_fe_index(range);
-        FECellIntegrator phi(matrix_free, 0, 0, 0, i, i);
+        FECellIntegrator phi(matrix_free, range);
 
         for (unsigned int cell = range.first; cell < range.second; ++cell)
           {
@@ -206,8 +205,7 @@ public:
   {
     matrix_free.template cell_loop<VectorType, VectorType>(
       [&](const auto &data, auto &dst, const auto &src, const auto range) {
-        const auto       i = data.get_cell_active_fe_index(range);
-        FECellIntegrator phi(matrix_free, 0, 0, 0, i, i);
+        FECellIntegrator phi(matrix_free, range);
 
         for (unsigned int cell = range.first; cell < range.second; ++cell)
           {
-- 
2.39.5