]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Simplify extract of face information, preparing for use in evaluation_kernels
authorMartin Kronbichler <martin.kronbichler@it.uu.se>
Fri, 10 Dec 2021 10:46:12 +0000 (11:46 +0100)
committerMartin Kronbichler <martin.kronbichler@it.uu.se>
Sat, 11 Dec 2021 11:14:21 +0000 (12:14 +0100)
examples/step-76/step-76.cc
include/deal.II/fe/mapping_q_internal.h
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/fe_evaluation_base_data.h
include/deal.II/matrix_free/mapping_info.templates.h
tests/matrix_free/interpolate_quadrature_01.cc

index c7b2a8f987ec7b1ee2c93f37fe0bf8e416487d00..dd437c18671a978d4c0c604222281c85af91ebb2 100644 (file)
@@ -790,11 +790,10 @@ namespace Euler_DG
                                                      VectorizedArrayType>::
                   template interpolate_quadrature<true, false>(
                     dim + 2,
+                    EvaluationFlags::values,
                     data.get_shape_info(),
                     buffer.data(),
                     phi_m.begin_values(),
-                    false,
-                    false,
                     face);
 
                 // Check if the face is an internal or a boundary face and
@@ -893,11 +892,10 @@ namespace Euler_DG
                                                      VectorizedArrayType>::
                   template interpolate_quadrature<false, true>(
                     dim + 2,
+                    EvaluationFlags::values,
                     data.get_shape_info(),
                     phi_m.begin_values(),
                     phi.begin_values(),
-                    false,
-                    false,
                     face);
               }
 
index 3d40e5154d7780e3c96dbe6088df4689a85cf3b6..1fed26230cda4f767d91dc3b217732635c5195a0 100644 (file)
@@ -1136,7 +1136,7 @@ namespace internal
       temp_data.descriptor.resize(1);
       temp_data.descriptor[0].n_q_points = n_q_points;
       FEEvaluationBaseData<dim, double, false, VectorizedArrayType> eval(
-        {&data.shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0);
+        {&data.shape_info, nullptr, &temp_data, 0, 0});
 
       // prepare arrays
       if (evaluation_flag != EvaluationFlags::nothing)
index 121f5a83d58f7c7b1c7e87631c768e2fb17ed65d..576187616bbb0db4124321f542f28fe6d40d9820 100644 (file)
@@ -2553,7 +2553,7 @@ namespace internal
 
   private:
     template <bool do_evaluate, bool add_into_output, int face_direction = 0>
-    static void
+    static DEAL_II_ALWAYS_INLINE void
     interpolate_generic(const unsigned int                     n_components,
                         const Number *                         input,
                         Number *                               output,
@@ -3137,6 +3137,20 @@ namespace internal
     // re-orientation
     std::array<const unsigned int *, n_face_orientations> orientation = {};
 
+#  ifdef DEBUG
+      // currently on structured meshes are supported -> face numbers and
+      // orientations have to be the same for all filled lanes
+      for (unsigned int v = 1; v < n_lanes; ++v)
+        {
+          if (this->all_face_numbers[v] != static_cast<std::uint8_t>(-1))
+            AssertDimension(this->all_face_numbers[0],
+                            this->all_face_numbers[v]);
+          if (this->all_face_orientations[v] != static_cast<std::uint8_t>(-1))
+            AssertDimension(this->all_face_orientations[0],
+                            this->all_face_orientations[v]);
+        }
+#  endif
+
     if (n_face_orientations == 1)
       orientation[0] =
         &eval.get_orientation_map()(eval.get_face_orientation(), 0);
index 81f6aaa3d4ee13591792b7162a2523f485611cf4..a3192566af9c74754d9f20439b84185b51b8d03a 100644 (file)
@@ -41,6 +41,8 @@
 #include <deal.II/matrix_free/type_traits.h>
 #include <deal.II/matrix_free/vector_access_internal.h>
 
+#include <type_traits>
+
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -615,26 +617,6 @@ public:
     AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
     const std::array<T, VectorizedArrayType::size()> &         value) const;
 
-  /**
-   * Return the id of the cells this FEEvaluation or FEFaceEvaluation is
-   * associated with.
-   */
-  std::array<unsigned int, VectorizedArrayType::size()>
-  get_cell_ids() const;
-
-  /**
-   * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is
-   * associated with.
-   */
-  std::array<unsigned int, VectorizedArrayType::size()>
-  get_cell_or_face_ids() const;
-
-  /**
-   * Return the first selected component.
-   */
-  unsigned int
-  get_first_selected_component() const;
-
   /**
    * Return the underlying MatrixFree object.
    */
@@ -792,18 +774,6 @@ protected:
    */
   const MatrixFree<dim, Number, VectorizedArrayType> *matrix_info;
 
-  /**
-   * Stores the number of components in the finite element as detected in the
-   * MatrixFree storage class for comparison with the template argument.
-   */
-  const unsigned int n_fe_components;
-
-  /**
-   * For a FiniteElement with more than one base element, select at which
-   * component this data structure should start.
-   */
-  const unsigned int first_selected_component;
-
   /**
    * A temporary data structure necessary to read degrees of freedom when no
    * MatrixFree object was given at initialization.
@@ -2802,20 +2772,6 @@ public:
    * points is inaccurate and this value must be used instead.
    */
   const unsigned int n_q_points;
-
-
-private:
-  /**
-   * Return face number of each face of the current face batch.
-   */
-  std::array<unsigned int, VectorizedArrayType::size()>
-  compute_face_no_data();
-
-  /**
-   * Determine the orientation of each face of the current face batch.
-   */
-  std::array<unsigned int, VectorizedArrayType::size()>
-  compute_face_orientations();
 };
 
 
@@ -2935,13 +2891,12 @@ inline FEEvaluationBase<dim,
                                                     n_q_points,
                                                     active_fe_index,
                                                     active_quad_index),
-      quad_no,
       is_interior_face,
-      face_type)
+      quad_no,
+      face_type,
+      first_selected_component)
   , scratch_data_array(data_in.acquire_scratch_data())
   , matrix_info(&data_in)
-  , n_fe_components(this->dof_info->start_components.back())
-  , first_selected_component(first_selected_component)
 {
   this->set_data_pointers(scratch_data_array, n_components_);
   Assert(
@@ -3002,13 +2957,11 @@ inline FEEvaluationBase<dim,
             MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
           mapping,
           quadrature,
-          update_flags))
+          update_flags),
+      n_components_,
+      first_selected_component)
   , scratch_data_array(new AlignedVector<VectorizedArrayType>())
   , matrix_info(nullptr)
-  , n_fe_components(n_components_)
-  // keep the number of the selected component within the current base element
-  // for reading dof values
-  , first_selected_component(first_selected_component)
 {
   const unsigned int base_element_number =
     fe.component_to_base_index(first_selected_component).first;
@@ -3052,8 +3005,6 @@ inline FEEvaluationBase<dim,
                          new AlignedVector<VectorizedArrayType>() :
                          other.matrix_info->acquire_scratch_data())
   , matrix_info(other.matrix_info)
-  , n_fe_components(other.n_fe_components)
-  , first_selected_component(other.first_selected_component)
 {
   if (other.matrix_info == nullptr)
     {
@@ -3143,9 +3094,6 @@ operator=(const FEEvaluationBase<dim,
       scratch_data_array = matrix_info->acquire_scratch_data();
     }
 
-  AssertDimension(n_fe_components, other.n_fe_components);
-  AssertDimension(first_selected_component, other.first_selected_component);
-
   this->set_data_pointers(scratch_data_array, n_components_);
 
   return *this;
@@ -3198,86 +3146,6 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 }
 
 
-
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-inline std::array<unsigned int, VectorizedArrayType::size()>
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  get_cell_ids() const
-{
-  Assert(this->matrix_info != nullptr, ExcNotInitialized());
-
-  const unsigned int                n_lanes = VectorizedArrayType::size();
-  std::array<unsigned int, n_lanes> cells;
-
-  // initialize array
-  for (unsigned int i = 0; i < n_lanes; ++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 < n_lanes; ++i)
-        cells[i] = this->cell * n_lanes + 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 < n_lanes; ++i)
-        {
-          // compute actual (non vectorized) cell ID
-          const unsigned int cell_this = this->cell * n_lanes + 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 / n_lanes)
-                          .cells_interior[face_index % n_lanes];
-          auto cell_p = this->matrix_info->get_face_info(face_index / n_lanes)
-                          .cells_exterior[face_index % n_lanes];
-
-          // 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_ =
-        this->is_interior_face ?
-          &this->matrix_info->get_face_info(this->cell).cells_interior[0] :
-          &this->matrix_info->get_face_info(this->cell).cells_exterior[0];
-      for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
-        if (cells_[i] != numbers::invalid_unsigned_int)
-          cells[i] = cells_[i];
-    }
-
-  return cells;
-}
-
-
 namespace internal
 {
   template <int dim,
@@ -3319,66 +3187,6 @@ namespace internal
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-std::array<unsigned int, VectorizedArrayType::size()>
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  get_cell_or_face_ids() const
-{
-  const unsigned int v_len = VectorizedArrayType::size();
-  std::array<unsigned int, VectorizedArrayType::size()> cells;
-
-  // initialize array
-  for (unsigned int i = 0; i < v_len; ++i)
-    cells[i] = numbers::invalid_unsigned_int;
-
-  if ((is_face &&
-       this->dof_access_index ==
-         internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-       this->is_interior_face == false) ||
-      (!is_face && !this->is_interior_face))
-    {
-      // cell-based face-loop: plus face
-      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 fn =
-            this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
-                                                                  this->face_no,
-                                                                  i);
-
-          if (fn == 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(fn / v_len)
-                          .cells_interior[fn % v_len];
-          auto cell_p = this->matrix_info->get_face_info(fn / v_len)
-                          .cells_exterior[fn % 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
-    {
-      for (unsigned int i = 0; i < v_len; ++i)
-        cells[i] = this->cell * v_len + i;
-    }
-
-  return cells;
-}
-
-
-
 template <int dim,
           int n_components_,
           typename Number,
@@ -3563,7 +3371,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
   Assert(this->dof_info != nullptr, ExcNotInitialized());
   Assert(this->matrix_info->indices_initialized() == true, ExcNotInitialized());
-  if (n_fe_components == 1)
+  if (this->n_fe_components == 1)
     for (unsigned int comp = 0; comp < n_components; ++comp)
       {
         Assert(src[comp] != nullptr,
@@ -3613,10 +3421,9 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     {
       for (unsigned int v = 0; v < n_vectorization_actual; ++v)
         if (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
-            this->dof_info
-                ->hanging_node_constraint_masks[(this->cell * n_lanes + v) *
-                                                  n_fe_components +
-                                                first_selected_component] !=
+            this->dof_info->hanging_node_constraint_masks
+                [(this->cell * n_lanes + v) * this->n_fe_components +
+                 this->first_selected_component] !=
               internal::MatrixFreeFunctions::ConstraintKinds::unconstrained)
           has_hn_constraints = true;
     }
@@ -3640,13 +3447,13 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     {
       const unsigned int *dof_indices =
         this->dof_info->dof_indices_interleaved.data() +
-        this->dof_info->row_starts[this->cell * n_fe_components * n_lanes]
+        this->dof_info->row_starts[this->cell * this->n_fe_components * n_lanes]
           .first +
         this->dof_info
             ->component_dof_indices_offset[this->active_fe_index]
                                           [this->first_selected_component] *
           n_lanes;
-      if (n_components == 1 || n_fe_components == 1)
+      if (n_components == 1 || this->n_fe_components == 1)
         for (unsigned int i = 0; i < dofs_per_component;
              ++i, dof_indices += n_lanes)
           for (unsigned int comp = 0; comp < n_components; ++comp)
@@ -3670,8 +3477,9 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   // to the dof indices of the cells on all lanes
   unsigned int        cells_copied[n_lanes];
   const unsigned int *cells;
-  bool                has_constraints  = false;
-  const unsigned int n_components_read = n_fe_components > 1 ? n_components : 1;
+  bool                has_constraints = false;
+  const unsigned int  n_components_read =
+    this->n_fe_components > 1 ? n_components : 1;
 
   if (is_face)
     {
@@ -3691,8 +3499,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
           Assert(cells[v] < this->dof_info->row_starts.size() - 1,
                  ExcInternalError());
           const std::pair<unsigned int, unsigned int> *my_index_start =
-            &this->dof_info->row_starts[cells[v] * n_fe_components +
-                                        first_selected_component];
+            &this->dof_info->row_starts[cells[v] * this->n_fe_components +
+                                        this->first_selected_component];
 
           // check whether any of the SIMD lanes has constraints, i.e., the
           // constraint indicator which is the second entry of row_starts
@@ -3709,23 +3517,22 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     }
   else
     {
-      AssertIndexRange((this->cell + 1) * n_lanes * n_fe_components,
+      AssertIndexRange((this->cell + 1) * n_lanes * this->n_fe_components,
                        this->dof_info->row_starts.size());
       for (unsigned int v = 0; v < n_vectorization_actual; ++v)
         {
           const std::pair<unsigned int, unsigned int> *my_index_start =
             &this->dof_info
-               ->row_starts[(this->cell * n_lanes + v) * n_fe_components +
-                            first_selected_component];
+               ->row_starts[(this->cell * n_lanes + v) * this->n_fe_components +
+                            this->first_selected_component];
           if (my_index_start[n_components_read].second !=
               my_index_start[0].second)
             has_constraints = true;
 
           if (this->dof_info->hanging_node_constraint_masks.size() > 0 &&
-              this->dof_info
-                  ->hanging_node_constraint_masks[(this->cell * n_lanes + v) *
-                                                    n_fe_components +
-                                                  first_selected_component] !=
+              this->dof_info->hanging_node_constraint_masks
+                  [(this->cell * n_lanes + v) * this->n_fe_components +
+                   this->first_selected_component] !=
                 internal::MatrixFreeFunctions::ConstraintKinds::unconstrained)
             has_hn_constraints = true;
 
@@ -3750,7 +3557,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
         for (unsigned int comp = 0; comp < n_components; ++comp)
           for (unsigned int i = 0; i < dofs_per_component; ++i)
             operation.process_empty(values_dofs[comp][i]);
-      if (n_components == 1 || n_fe_components == 1)
+      if (n_components == 1 || this->n_fe_components == 1)
         {
           for (unsigned int v = 0; v < n_vectorization_actual; ++v)
             for (unsigned int i = 0; i < dofs_per_component; ++i)
@@ -3786,9 +3593,9 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
       const unsigned int cell_index =
         is_face ? cells[v] : this->cell * n_lanes + v;
       const unsigned int cell_dof_index =
-        cell_index * n_fe_components + first_selected_component;
+        cell_index * this->n_fe_components + this->first_selected_component;
       const unsigned int n_components_read =
-        n_fe_components > 1 ? n_components : 1;
+        this->n_fe_components > 1 ? n_components : 1;
       unsigned int index_indicators =
         this->dof_info->row_starts[cell_dof_index].second;
       unsigned int next_index_indicators =
@@ -3816,7 +3623,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
           next_index_indicators = index_indicators;
         }
 
-      if (n_components == 1 || n_fe_components == 1)
+      if (n_components == 1 || this->n_fe_components == 1)
         {
           unsigned int ind_local = 0;
           for (; index_indicators != next_index_indicators; ++index_indicators)
@@ -3950,7 +3757,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   Assert(!local_dof_indices.empty(), ExcNotInitialized());
 
   const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
-  unsigned int      index = first_selected_component * dofs_per_component;
+  unsigned int      index = this->first_selected_component * dofs_per_component;
   for (unsigned int comp = 0; comp < n_components; ++comp)
     {
       for (unsigned int i = 0; i < dofs_per_component; ++i, ++index)
@@ -4021,10 +3828,11 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     {
       const unsigned int dof_index =
         dof_indices_cont[this->cell * VectorizedArrayType::size()] +
-        this->dof_info->component_dof_indices_offset[this->active_fe_index]
-                                                    [first_selected_component] *
+        this->dof_info
+            ->component_dof_indices_offset[this->active_fe_index]
+                                          [this->first_selected_component] *
           VectorizedArrayType::size();
-      if (n_components == 1 || n_fe_components == 1)
+      if (n_components == 1 || this->n_fe_components == 1)
         for (unsigned int comp = 0; comp < n_components; ++comp)
           operation.process_dofs_vectorized(dofs_per_component,
                                             dof_index,
@@ -4040,7 +3848,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
       return;
     }
 
-  std::array<unsigned int, VectorizedArrayType::size()> cells =
+  const std::array<unsigned int, VectorizedArrayType::size()> &cells =
     this->get_cell_or_face_ids();
 
   // More general case: Must go through the components one by one and apply
@@ -4102,7 +3910,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
       if (n_filled_lanes == VectorizedArrayType::size() &&
           n_lanes == VectorizedArrayType::size() && !is_ecl)
         {
-          if (n_components == 1 || n_fe_components == 1)
+          if (n_components == 1 || this->n_fe_components == 1)
             {
               for (unsigned int comp = 0; comp < n_components; ++comp)
                 {
@@ -4128,12 +3936,12 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
         for (unsigned int comp = 0; comp < n_components; ++comp)
           {
             auto vector_ptrs = compute_vector_ptrs(
-              (n_components == 1 || n_fe_components == 1) ? comp : 0);
+              (n_components == 1 || this->n_fe_components == 1) ? comp : 0);
 
             for (unsigned int i = 0; i < dofs_per_component; ++i)
               operation.process_empty(values_dofs[comp][i]);
 
-            if (n_components == 1 || n_fe_components == 1)
+            if (n_components == 1 || this->n_fe_components == 1)
               {
                 for (unsigned int v = 0; v < n_filled_lanes; ++v)
                   if (mask[v] == true)
@@ -4179,7 +3987,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
           internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
             contiguous)
         {
-          if (n_components == 1 || n_fe_components == 1)
+          if (n_components == 1 || this->n_fe_components == 1)
             for (unsigned int comp = 0; comp < n_components; ++comp)
               operation.process_dofs_vectorized_transpose(dofs_per_component,
                                                           dof_indices,
@@ -4198,7 +4006,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
                  interleaved_contiguous_strided)
         {
-          if (n_components == 1 || n_fe_components == 1)
+          if (n_components == 1 || this->n_fe_components == 1)
             for (unsigned int i = 0; i < dofs_per_component; ++i)
               {
                 for (unsigned int comp = 0; comp < n_components; ++comp)
@@ -4229,7 +4037,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
           const unsigned int *offsets =
             &this->dof_info->dof_indices_interleave_strides
                [ind][VectorizedArrayType::size() * this->cell];
-          if (n_components == 1 || n_fe_components == 1)
+          if (n_components == 1 || this->n_fe_components == 1)
             for (unsigned int i = 0; i < dofs_per_component; ++i)
               {
                 for (unsigned int comp = 0; comp < n_components; ++comp)
@@ -4266,7 +4074,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
               contiguous)
           {
-            if (n_components == 1 || n_fe_components == 1)
+            if (n_components == 1 || this->n_fe_components == 1)
               {
                 for (unsigned int v = 0; v < n_filled_lanes; ++v)
                   if (mask[v] == true)
@@ -4293,7 +4101,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                  [ind][VectorizedArrayType::size() * this->cell];
             for (unsigned int v = 0; v < n_filled_lanes; ++v)
               AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
-            if (n_components == 1 || n_fe_components == 1)
+            if (n_components == 1 || this->n_fe_components == 1)
               for (unsigned int v = 0; v < n_filled_lanes; ++v)
                 {
                   if (mask[v] == true)
@@ -4472,7 +4280,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
       const unsigned int cell_index =
         is_face ? cells[v] : this->cell * n_lanes + v;
       const unsigned int cell_dof_index =
-        cell_index * n_fe_components + first_selected_component;
+        cell_index * this->n_fe_components + this->first_selected_component;
 
       const auto mask =
         this->dof_info->hanging_node_constraint_masks[cell_dof_index];
@@ -4662,18 +4470,6 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
 /*------------------------------ access to data fields ----------------------*/
 
-template <int dim,
-          int n_components,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-inline unsigned int
-FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
-  get_first_selected_component() const
-{
-  return first_selected_component;
-}
-
 
 
 template <int dim,
@@ -7062,6 +6858,9 @@ FEEvaluation<dim,
   this->jacobian = &this->mapping_data->jacobians[0][offsets];
   this->J_value  = &this->mapping_data->JxW_values[offsets];
 
+  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
+    this->cell_ids[i] = cell_index * VectorizedArrayType::size() + i;
+
 #  ifdef DEBUG
   this->dof_values_initialized     = false;
   this->values_quad_initialized    = false;
@@ -7353,7 +7152,6 @@ namespace internal
    * Implementation for standard vectors (that have the begin() methods).
    */
   template <typename Number,
-            typename VectorizedArrayType,
             typename VectorType,
             typename T,
             typename std::enable_if<
@@ -7363,26 +7161,27 @@ namespace internal
               VectorType>::type * = nullptr>
   bool
   try_gather_evaluate_inplace(
-    T                                             phi,
-    const VectorType &                            input_vector,
-    const unsigned int                            cell,
-    const unsigned int                            active_fe_index,
-    const unsigned int                            first_selected_component,
-    const internal::MatrixFreeFunctions::DoFInfo *dof_info,
-    const EvaluationFlags::EvaluationFlags        evaluation_flag)
+    T &                                    phi,
+    const VectorType &                     input_vector,
+    const EvaluationFlags::EvaluationFlags evaluation_flag)
   {
+    const unsigned int cell     = phi.get_current_cell_index();
+    const auto &       dof_info = phi.get_dof_info();
+    using VectorizedArrayType =
+      typename std::remove_reference<decltype(*phi.begin_dof_values())>::type;
+
     // If the index storage is interleaved and contiguous and the vector storage
     // has the correct alignment, we can directly pass the pointer into the
     // vector to the evaluate() call, without reading the vector entries into a
     // separate data field. This saves some operations.
     if (std::is_same<typename VectorType::value_type, Number>::value &&
-        dof_info->index_storage_variants
+        dof_info.index_storage_variants
             [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] ==
           internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
             interleaved_contiguous &&
         reinterpret_cast<std::size_t>(
           input_vector.begin() +
-          dof_info->dof_indices_contiguous
+          dof_info.dof_indices_contiguous
             [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
             [cell * VectorizedArrayType::size()]) %
             sizeof(VectorizedArrayType) ==
@@ -7391,14 +7190,15 @@ namespace internal
         const VectorizedArrayType *vec_values =
           reinterpret_cast<const VectorizedArrayType *>(
             input_vector.begin() +
-            dof_info->dof_indices_contiguous
+            dof_info.dof_indices_contiguous
               [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
               [cell * VectorizedArrayType::size()] +
-            dof_info->component_dof_indices_offset[active_fe_index]
-                                                  [first_selected_component] *
+            dof_info.component_dof_indices_offset
+                [phi.get_active_fe_index()]
+                [phi.get_first_selected_component()] *
               VectorizedArrayType::size());
 
-        phi->evaluate(vec_values, evaluation_flag);
+        phi.evaluate(vec_values, evaluation_flag);
 
         return true;
       }
@@ -7410,7 +7210,6 @@ namespace internal
    * Implementation for block vectors.
    */
   template <typename Number,
-            typename VectorizedArrayType,
             typename VectorType,
             typename T,
             typename std::enable_if<
@@ -7419,12 +7218,8 @@ namespace internal
                               Number *>::value,
               VectorType>::type * = nullptr>
   bool
-  try_gather_evaluate_inplace(T,
+  try_gather_evaluate_inplace(T &,
                               const VectorType &,
-                              const unsigned int,
-                              const unsigned int,
-                              const unsigned int,
-                              const internal::MatrixFreeFunctions::DoFInfo *,
                               const EvaluationFlags::EvaluationFlags)
   {
     return false;
@@ -7434,7 +7229,6 @@ namespace internal
    * Implementation for vectors that have the begin() methods.
    */
   template <typename Number,
-            typename VectorizedArrayType,
             typename VectorType,
             typename T,
             typename std::enable_if<
@@ -7444,26 +7238,27 @@ namespace internal
               VectorType>::type * = nullptr>
   bool
   try_integrate_scatter_inplace(
-    T                                             phi,
-    VectorType &                                  destination,
-    const unsigned int                            cell,
-    const unsigned int                            active_fe_index,
-    const unsigned int                            first_selected_component,
-    const internal::MatrixFreeFunctions::DoFInfo *dof_info,
-    const EvaluationFlags::EvaluationFlags        evaluation_flag)
+    T &                                    phi,
+    VectorType &                           destination,
+    const EvaluationFlags::EvaluationFlags evaluation_flag)
   {
+    const unsigned int cell     = phi.get_current_cell_index();
+    const auto &       dof_info = phi.get_dof_info();
+    using VectorizedArrayType =
+      typename std::remove_reference<decltype(*phi.begin_dof_values())>::type;
+
     // If the index storage is interleaved and contiguous and the vector storage
     // has the correct alignment, we can directly pass the pointer into the
     // vector to the evaluate() call, without reading the vector entries into a
     // separate data field. This saves some operations.
     if (std::is_same<typename VectorType::value_type, Number>::value &&
-        dof_info->index_storage_variants
+        dof_info.index_storage_variants
             [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell][cell] ==
           internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
             interleaved_contiguous &&
         reinterpret_cast<std::size_t>(
           destination.begin() +
-          dof_info->dof_indices_contiguous
+          dof_info.dof_indices_contiguous
             [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
             [cell * VectorizedArrayType::size()]) %
             sizeof(VectorizedArrayType) ==
@@ -7475,11 +7270,12 @@ namespace internal
             dof_info->dof_indices_contiguous
               [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
               [cell * VectorizedArrayType::size()] +
-            dof_info->component_dof_indices_offset[active_fe_index]
-                                                  [first_selected_component] *
+            dof_info->component_dof_indices_offset
+                [phi.get_active_fe_index()]
+                [phi.get_first_selected_component()] *
               VectorizedArrayType::size());
 
-        phi->integrate(evaluation_flag, vec_values, true);
+        phi.integrate(evaluation_flag, vec_values, true);
 
         return true;
       }
@@ -7491,7 +7287,6 @@ namespace internal
    * Implementation for all other vectors like block vectors.
    */
   template <typename Number,
-            typename VectorizedArrayType,
             typename VectorType,
             typename T,
             typename std::enable_if<
@@ -7502,10 +7297,6 @@ namespace internal
   bool
   try_integrate_scatter_inplace(T,
                                 VectorType &,
-                                const unsigned int,
-                                const unsigned int,
-                                const unsigned int,
-                                const internal::MatrixFreeFunctions::DoFInfo *,
                                 const EvaluationFlags::EvaluationFlags)
   {
     return false;
@@ -7531,14 +7322,9 @@ FEEvaluation<dim,
   gather_evaluate(const VectorType &                     input_vector,
                   const EvaluationFlags::EvaluationFlags evaluation_flag)
 {
-  if (internal::try_gather_evaluate_inplace<Number, VectorizedArrayType>(
-        this,
-        input_vector,
-        this->cell,
-        this->active_fe_index,
-        this->first_selected_component,
-        this->dof_info,
-        evaluation_flag) == false)
+  if (internal::try_gather_evaluate_inplace<Number>(*this,
+                                                    input_vector,
+                                                    evaluation_flag) == false)
     {
       this->read_dof_values(input_vector);
       evaluate(this->begin_dof_values(), evaluation_flag);
@@ -7743,14 +7529,8 @@ FEEvaluation<dim,
   integrate_scatter(const EvaluationFlags::EvaluationFlags integration_flag,
                     VectorType &                           destination)
 {
-  if (internal::try_integrate_scatter_inplace<Number, VectorizedArrayType>(
-        this,
-        destination,
-        this->cell,
-        this->active_fe_index,
-        this->first_selected_component,
-        this->dof_info,
-        integration_flag) == false)
+  if (internal::try_integrate_scatter_inplace<Number>(
+        *this, destination, integration_flag) == false)
     {
       integrate(integration_flag, this->begin_dof_values());
       this->distribute_local_to_global(destination);
@@ -7871,6 +7651,8 @@ FEFaceEvaluation<dim,
              "is_interior_face set to true. "));
 
   this->reinit_face(this->matrix_info->get_face_info(face_index));
+  for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
+    this->cell_or_face_ids[i] = face_index * VectorizedArrayType::size() + i;
 
   this->cell_type = this->matrix_info->get_mapping_info().face_type[face_index];
   const unsigned int offsets =
@@ -7925,12 +7707,80 @@ FEFaceEvaluation<dim,
 
   this->cell_type = this->matrix_info->get_mapping_info().cell_type[cell_index];
   this->cell      = cell_index;
-  this->face_orientation = 0;
-  this->subface_index    = GeometryInfo<dim>::max_children_per_cell;
-  this->face_no          = face_number;
+  this->subface_index = GeometryInfo<dim>::max_children_per_cell;
   this->dof_access_index =
     internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
 
+  constexpr unsigned int n_lanes = VectorizedArrayType::size();
+
+  if (this->is_interior_face == false)
+    {
+      // 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 < n_lanes; ++i)
+        {
+          // compute actual (non vectorized) cell ID
+          const unsigned int cell_this = cell_index * n_lanes + i;
+          // compute face ID
+          unsigned int face_index =
+            this->matrix_info->get_cell_and_face_to_plain_faces()(cell_index,
+                                                                  face_number,
+                                                                  i);
+
+          if (face_index == numbers::invalid_unsigned_int)
+            {
+              this->cell_ids[i]              = numbers::invalid_unsigned_int;
+              this->all_face_numbers[i]      = numbers::invalid_unsigned_int;
+              this->all_face_orientations[i] = numbers::invalid_unsigned_int;
+              continue; // invalid face ID: no neighbor on boundary
+            }
+
+          const auto &faces =
+            this->matrix_info->get_face_info(face_index / n_lanes);
+          // get cell ID on both sides of face
+          auto cell_m = faces.cells_interior[face_index % n_lanes];
+          auto cell_p = faces.cells_exterior[face_index % n_lanes];
+
+          const bool is_interior_face = cell_m != cell_this;
+
+          Assert(cell_m == cell_this || cell_p == cell_this,
+                 ExcInternalError());
+
+          // compare the IDs with the given cell ID
+          if (is_interior_face)
+            {
+              this->cell_ids[i]         = cell_m; // neighbor has the other ID
+              this->all_face_numbers[i] = faces.interior_face_no;
+            }
+          else
+            {
+              this->cell_ids[i]         = cell_p;
+              this->all_face_numbers[i] = faces.exterior_face_no;
+            }
+
+          const bool   orientation_interior_face = faces.face_orientation >= 8;
+          unsigned int face_orientation          = faces.face_orientation % 8;
+          if (is_interior_face != orientation_interior_face)
+            {
+              constexpr std::array<std::uint8_t, 8> table{
+                {0, 1, 0, 3, 6, 5, 4, 7}};
+              face_orientation = table[face_orientation];
+            }
+          this->all_face_orientations[i] = face_orientation;
+        }
+
+      this->face_no          = this->all_face_numbers[0];
+      this->face_orientation = this->all_face_orientations[0];
+    }
+  else
+    {
+      this->face_orientation = 0;
+      this->face_no          = face_number;
+      for (unsigned int i = 0; i < n_lanes; ++i)
+        this->cell_ids[i] = cell_index * n_lanes + i;
+    }
+
   const unsigned int offsets =
     this->matrix_info->get_mapping_info()
       .face_data_by_cells[this->quad_no]
@@ -8470,173 +8320,6 @@ FEFaceEvaluation<dim,
 
 
 
-template <int dim,
-          int fe_degree,
-          int n_q_points_1d,
-          int n_components_,
-          typename Number,
-          typename VectorizedArrayType>
-std::array<unsigned int, VectorizedArrayType::size()>
-FEFaceEvaluation<dim,
-                 fe_degree,
-                 n_q_points_1d,
-                 n_components_,
-                 Number,
-                 VectorizedArrayType>::compute_face_no_data()
-{
-  std::array<unsigned int, VectorizedArrayType::size()> face_no_data;
-
-  if (this->dof_access_index !=
-        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ||
-      this->is_interior_face == true)
-    {
-      std::fill(face_no_data.begin(),
-                face_no_data.begin() +
-                  this->dof_info->n_vectorization_lanes_filled
-                    [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-                    [this->cell],
-                this->face_no);
-    }
-  else
-    {
-      std::fill(face_no_data.begin(),
-                face_no_data.end(),
-                numbers::invalid_unsigned_int);
-
-      for (unsigned int i = 0;
-           i < this->dof_info->n_vectorization_lanes_filled
-                 [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-                 [this->cell];
-           ++i)
-        {
-          // compute actual (non vectorized) cell ID
-          const unsigned int cell_this =
-            this->cell * VectorizedArrayType::size() + i;
-          // compute face ID
-          const unsigned int face_index =
-            this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
-                                                                  this->face_no,
-                                                                  i);
-
-          Assert(face_index != numbers::invalid_unsigned_int,
-                 ExcNotInitialized());
-
-          // get cell ID on both sides of face
-          auto cell_m =
-            this->matrix_info
-              ->get_face_info(face_index / VectorizedArrayType::size())
-              .cells_interior[face_index % VectorizedArrayType::size()];
-
-          // compare the IDs with the given cell ID
-          face_no_data[i] =
-            (cell_m == cell_this) ?
-              this->matrix_info
-                ->get_face_info(face_index / VectorizedArrayType::size())
-                .exterior_face_no :
-              this->matrix_info
-                ->get_face_info(face_index / VectorizedArrayType::size())
-                .interior_face_no;
-        }
-    }
-
-  return face_no_data;
-}
-
-
-
-template <int dim,
-          int fe_degree,
-          int n_q_points_1d,
-          int n_components_,
-          typename Number,
-          typename VectorizedArrayType>
-std::array<unsigned int, VectorizedArrayType::size()>
-FEFaceEvaluation<dim,
-                 fe_degree,
-                 n_q_points_1d,
-                 n_components_,
-                 Number,
-                 VectorizedArrayType>::compute_face_orientations()
-{
-  std::array<unsigned int, VectorizedArrayType::size()> face_no_data;
-
-  if (this->dof_access_index !=
-        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ||
-      this->is_interior_face == true)
-    {
-      Assert(false, ExcNotImplemented());
-    }
-  else
-    {
-      std::fill(face_no_data.begin(),
-                face_no_data.end(),
-                numbers::invalid_unsigned_int);
-
-      if (dim == 3)
-        {
-          for (unsigned int i = 0;
-               i < this->dof_info->n_vectorization_lanes_filled
-                     [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-                     [this->cell];
-               ++i)
-            {
-              // compute actual (non vectorized) cell ID
-              const unsigned int cell_this =
-                this->cell * VectorizedArrayType::size() + i;
-              // compute face ID
-              const unsigned int face_index =
-                this->matrix_info->get_cell_and_face_to_plain_faces()(
-                  this->cell, this->face_no, i);
-
-              Assert(face_index != numbers::invalid_unsigned_int,
-                     ExcNotInitialized());
-
-              const unsigned int macro =
-                face_index / VectorizedArrayType::size();
-              const unsigned int lane =
-                face_index % VectorizedArrayType::size();
-
-              const auto &faces = this->matrix_info->get_face_info(macro);
-
-              // get cell ID on both sides of face
-              auto cell_m = faces.cells_interior[lane];
-
-              const bool is_interior_face = cell_m != cell_this;
-              const bool fo_interior_face = faces.face_orientation >= 8;
-
-              unsigned int face_orientation = faces.face_orientation % 8;
-
-              if (is_interior_face != fo_interior_face)
-                {
-                  // invert (see also:
-                  // Triangulation::update_periodic_face_map())
-                  static const std::array<unsigned int, 8> table{
-                    {0, 1, 2, 3, 6, 5, 4, 7}};
-
-                  face_orientation = table[face_orientation];
-                }
-
-              // compare the IDs with the given cell ID
-              face_no_data[i] = face_orientation;
-            }
-        }
-      else
-        {
-          std::fill(
-            face_no_data.begin(),
-            face_no_data.begin() +
-              this->dof_info->n_vectorization_lanes_filled
-                [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-                [this->cell],
-            0);
-        }
-    }
-
-  return face_no_data;
-}
-
-
-
 template <int dim,
           int fe_degree,
           int n_q_points_1d,
index bf7457a1e528db44d773cf429ed6f52f66d65e2f..447297c1c8a39b19f57e06ab84153be1bd402bc5 100644 (file)
@@ -84,8 +84,10 @@ namespace internal
  * This base class of the FEEvaluation and FEFaceEvaluation classes handles
  * mapping-related information independent of the degrees of freedom and
  * finite element in use. This class provides access functionality for user
- * code but is otherwise invisible without any public constructor. The usage
- * is through the class FEEvaluation instead.
+ * code. The main usage is through the class FEEvaluation. However, there is
+ * functionality to construct an object of type FEEvaluationBaseData given
+ * suitable data pointers to internal data, which allows usage in some
+ * scenarios where no full MatrixFree object is available.
  *
  * This class has four template arguments:
  *
@@ -137,9 +139,10 @@ public:
                                         const MappingInfoStorageType *,
                                         unsigned int,
                                         unsigned int> &info,
-                       const unsigned int              quad_no,
-                       const bool                      is_interior_face,
-                       const unsigned int              face_type);
+                       const bool                      is_interior_face = true,
+                       const unsigned int              quad_no          = 0,
+                       const unsigned int              face_type        = 0,
+                       const unsigned int first_selected_component      = 0);
 
   /**
    * Copy constructor.
@@ -404,6 +407,12 @@ public:
   unsigned int
   get_active_quadrature_index() const;
 
+  /**
+   * Return the first selected component in a multi-component system.
+   */
+  unsigned int
+  get_first_selected_component() const;
+
   /**
    * If is_face is true, this function returns the face number within a cell for
    * the face this object was initialized to. On cells where the face makes no
@@ -471,6 +480,46 @@ public:
   bool
   get_is_interior_face() const;
 
+  /**
+   * Return the id of the cells this FEEvaluation or FEFaceEvaluation is
+   * associated with.
+   */
+  const std::array<unsigned int, VectorizedArrayType::size()> &
+  get_cell_ids() const;
+
+  /**
+   * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is
+   * associated with.
+   */
+  const std::array<unsigned int, VectorizedArrayType::size()> &
+  get_cell_or_face_ids() const;
+
+  /**
+   * Return the (non-vectorized) number of faces within cells in case of ECL
+   * for the exterior cells where the single number accessible via
+   * get_face_no() is not enough to cover all cases.
+   *
+   * @note Only available for `dof_access_index == dof_access_cell` and
+   * `is_interior_face == false`.
+   *
+   * @note This function depends on the internal representation of data, which
+   * is not stable between releases of deal.II, and is hence mostly for
+   * internal use.
+   */
+  const std::array<std::uint8_t, VectorizedArrayType::size()> &
+  get_all_face_numbers() const;
+
+  /**
+   * Store the orientation of the neighbor's faces with respect to the current
+   * cell for the case of exterior faces on ECL with possibly different
+   * orientations behind different cells.
+   *
+   * @note Only available for `dof_access_index == dof_access_cell` and
+   * `is_interior_face == false`.
+   */
+  const std::array<std::uint8_t, VectorizedArrayType::size()> &
+  get_all_face_orientations() const;
+
   //@}
 
 protected:
@@ -481,7 +530,9 @@ protected:
   FEEvaluationBaseData(
     const std::shared_ptr<
       internal::MatrixFreeFunctions::
-        MappingDataOnTheFly<dim, Number, VectorizedArrayType>> &mapping_data);
+        MappingDataOnTheFly<dim, Number, VectorizedArrayType>> &mapping_data,
+    const unsigned int                                          n_fe_components,
+    const unsigned int first_selected_component);
 
   /**
    * A pointer to the unit cell shape data, i.e., values, gradients and
@@ -513,6 +564,18 @@ protected:
    */
   const unsigned int quad_no;
 
+  /**
+   * Stores the number of components in the finite element as detected in the
+   * MatrixFree storage class for comparison with the template argument.
+   */
+  const unsigned int n_fe_components;
+
+  /**
+   * For a FiniteElement with more than one base element, select at which
+   * component this data structure should start.
+   */
+  const unsigned int first_selected_component;
+
   /**
    * The active FE index for this class for efficient indexing in the hp-case.
    */
@@ -731,6 +794,38 @@ protected:
    */
   internal::MatrixFreeFunctions::GeometryType cell_type;
 
+  /**
+   * Stores the (non-vectorized) id of the cells or faces this object is
+   * initialized to. Relevant for ECL.
+   */
+  std::array<unsigned int, VectorizedArrayType::size()> cell_ids;
+
+  /**
+   * Stores the (non-vectorized) id of the cells or faces this object is
+   * initialized to. Relevant for ECL.
+   */
+  std::array<unsigned int, VectorizedArrayType::size()> cell_or_face_ids;
+
+  /**
+   * Stores the (non-vectorized) number of faces within cells in case of ECL
+   * for the exterior cells where the single number `face_no` is not enough to
+   * cover all cases.
+   *
+   * @note Only available for `dof_access_index == dof_access_cell` and
+   * `is_interior_face == false`.
+   */
+  std::array<std::uint8_t, VectorizedArrayType::size()> all_face_numbers;
+
+  /**
+   * Store the orientation of the neighbor's faces with respect to the current
+   * cell for the case of exterior faces on ECL with possibly different
+   * orientations behind different cells.
+   *
+   * @note Only available for `dof_access_index == dof_access_cell` and
+   * `is_interior_face == false`.
+   */
+  std::array<std::uint8_t, VectorizedArrayType::size()> all_face_orientations;
+
   /**
    * Geometry data that can be generated FEValues on the fly with the
    * respective constructor, as an alternative to the entry point with
@@ -762,13 +857,16 @@ inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
                                         const MappingInfoStorageType *,
                                         unsigned int,
                                         unsigned int> &shape_dof_info,
-                       const unsigned int              quad_no,
                        const bool                      is_interior_face,
-                       const unsigned int              face_type)
+                       const unsigned int              quad_no,
+                       const unsigned int              face_type,
+                       const unsigned int              first_selected_component)
   : data(std::get<0>(shape_dof_info))
   , dof_info(std::get<1>(shape_dof_info))
   , mapping_data(std::get<2>(shape_dof_info))
   , quad_no(quad_no)
+  , n_fe_components(dof_info != nullptr ? dof_info->start_components.back() : 0)
+  , first_selected_component(first_selected_component)
   , active_fe_index(std::get<3>(shape_dof_info))
   , active_quad_index(std::get<4>(shape_dof_info))
   , descriptor(
@@ -810,11 +908,15 @@ inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
   FEEvaluationBaseData(
     const std::shared_ptr<
       internal::MatrixFreeFunctions::
-        MappingDataOnTheFly<dim, Number, VectorizedArrayType>> &mapped_geometry)
+        MappingDataOnTheFly<dim, Number, VectorizedArrayType>> &mapped_geometry,
+    const unsigned int                                          n_fe_components,
+    const unsigned int first_selected_component)
   : data(nullptr)
   , dof_info(nullptr)
   , mapping_data(nullptr)
   , quad_no(numbers::invalid_unsigned_int)
+  , n_fe_components(n_fe_components)
+  , first_selected_component(first_selected_component)
   , active_fe_index(numbers::invalid_unsigned_int)
   , active_quad_index(numbers::invalid_unsigned_int)
   , descriptor(nullptr)
@@ -844,6 +946,8 @@ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::operator=(
   const FEEvaluationBaseData &other)
 {
   AssertDimension(quad_no, other.quad_no);
+  AssertDimension(n_fe_components, other.n_fe_components);
+  AssertDimension(first_selected_component, other.first_selected_component);
   AssertDimension(active_fe_index, other.active_fe_index);
   AssertDimension(active_quad_index, other.active_quad_index);
   AssertDimension(n_quadrature_points, descriptor->n_q_points);
@@ -883,6 +987,78 @@ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::operator=(
 
 
 
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline void
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  set_data_pointers(AlignedVector<VectorizedArrayType> *scratch_data_array,
+                    const unsigned int                  n_components)
+{
+  Assert(scratch_data_array != nullptr, ExcInternalError());
+
+  const unsigned int tensor_dofs_per_component =
+    Utilities::fixed_power<dim>(data->data.front().fe_degree + 1);
+  const unsigned int dofs_per_component = data->dofs_per_component_on_cell;
+
+  const unsigned int size_scratch_data =
+    std::max(tensor_dofs_per_component + 1, dofs_per_component) * n_components *
+      3 +
+    2 * n_quadrature_points;
+  const unsigned int size_data_arrays =
+    n_components * dofs_per_component +
+    (n_components * ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
+     n_quadrature_points);
+
+  const unsigned int allocated_size = size_scratch_data + size_data_arrays;
+  scratch_data_array->resize_fast(allocated_size);
+  scratch_data.reinit(scratch_data_array->begin() + size_data_arrays,
+                      size_scratch_data);
+
+  // set the pointers to the correct position in the data array
+  values_dofs = scratch_data_array->begin();
+  values_quad = scratch_data_array->begin() + n_components * dofs_per_component;
+  gradients_quad = scratch_data_array->begin() +
+                   n_components * (dofs_per_component + n_quadrature_points);
+  gradients_from_hessians_quad =
+    scratch_data_array->begin() +
+    n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
+  hessians_quad =
+    scratch_data_array->begin() +
+    n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points);
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline void
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::reinit_face(
+  const internal::MatrixFreeFunctions::FaceToCellTopology<
+    VectorizedArrayType::size()> &face)
+{
+  Assert(is_face == true,
+         ExcMessage("Faces can only be set if the is_face template parameter "
+                    "is true"));
+  face_no = (is_interior_face ? face.interior_face_no : face.exterior_face_no);
+  subface_index = is_interior_face == true ?
+                    GeometryInfo<dim>::max_children_per_cell :
+                    face.subface_index;
+
+  // First check if interior or exterior cell has non-standard orientation
+  // (i.e. the third bit is one or not). Then set zero if this cell has
+  // standard-orientation else copy the first three bits
+  // (which is equivalent to modulo 8). See also the documentation of
+  // internal::MatrixFreeFunctions::FaceToCellTopology::face_orientation.
+  face_orientation = (is_interior_face == (face.face_orientation >= 8)) ?
+                       (face.face_orientation % 8) :
+                       0;
+
+  if (is_interior_face)
+    cell_ids = face.cells_interior;
+  else
+    cell_ids = face.cells_exterior;
+}
+
+
+
 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
 FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
@@ -1147,6 +1323,16 @@ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
 
 
 
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_first_selected_component() const
+{
+  return first_selected_component;
+}
+
+
+
 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline ArrayView<VectorizedArrayType>
 FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
@@ -1219,65 +1405,63 @@ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
 
 
 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline void
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  set_data_pointers(AlignedVector<VectorizedArrayType> *scratch_data_array,
-                    const unsigned int                  n_components)
+inline const std::array<unsigned int, VectorizedArrayType::size()> &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_ids()
+  const
 {
-  Assert(scratch_data_array != nullptr, ExcInternalError());
-
-  const unsigned int tensor_dofs_per_component =
-    Utilities::fixed_power<dim>(data->data.front().fe_degree + 1);
-  const unsigned int dofs_per_component = data->dofs_per_component_on_cell;
+  Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  return cell_ids;
+}
 
-  const unsigned int size_scratch_data =
-    std::max(tensor_dofs_per_component + 1, dofs_per_component) * n_components *
-      3 +
-    2 * n_quadrature_points;
-  const unsigned int size_data_arrays =
-    n_components * dofs_per_component +
-    (n_components * ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
-     n_quadrature_points);
 
-  const unsigned int allocated_size = size_scratch_data + size_data_arrays;
-  scratch_data_array->resize_fast(allocated_size);
-  scratch_data.reinit(scratch_data_array->begin() + size_data_arrays,
-                      size_scratch_data);
 
-  // set the pointers to the correct position in the data array
-  values_dofs = scratch_data_array->begin();
-  values_quad = scratch_data_array->begin() + n_components * dofs_per_component;
-  gradients_quad = scratch_data_array->begin() +
-                   n_components * (dofs_per_component + n_quadrature_points);
-  gradients_from_hessians_quad =
-    scratch_data_array->begin() +
-    n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
-  hessians_quad =
-    scratch_data_array->begin() +
-    n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points);
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline const std::array<unsigned int, VectorizedArrayType::size()> &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_cell_or_face_ids() const
+{
+  Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  if (!is_face || dof_access_index ==
+                    internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
+    return cell_ids;
+  else
+    return cell_or_face_ids;
 }
 
 
 
 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline void
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::reinit_face(
-  const internal::MatrixFreeFunctions::FaceToCellTopology<
-    VectorizedArrayType::size()> &face)
+inline const std::array<std::uint8_t, VectorizedArrayType::size()> &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_all_face_numbers() const
 {
-  face_no = (is_interior_face ? face.interior_face_no : face.exterior_face_no);
-  subface_index = is_interior_face == true ?
-                    GeometryInfo<dim>::max_children_per_cell :
-                    face.subface_index;
+  Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(is_face &&
+           dof_access_index ==
+             internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+           is_interior_face == false,
+         ExcMessage("All face numbers can only be queried for ECL at exterior "
+                    "faces. Use get_face_no() in other cases."));
+
+  return all_face_numbers;
+}
 
-  // First check if interior or exterior cell has non-standard orientation
-  // (i.e. the third bit is one or not). Then set zero if this cell has
-  // standard-orientation else copy the first three bits
-  // (which is equivalent to modulo 8). See also the documentation of
-  // internal::MatrixFreeFunctions::FaceToCellTopology::face_orientation.
-  face_orientation = (is_interior_face == (face.face_orientation >= 8)) ?
-                       (face.face_orientation % 8) :
-                       0;
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline const std::array<std::uint8_t, VectorizedArrayType::size()> &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_all_face_orientations() const
+{
+  Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  Assert(is_face &&
+           dof_access_index ==
+             internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+           is_interior_face == false,
+         ExcMessage("All face numbers can only be queried for ECL at exterior "
+                    "faces. Use get_face_no() in other cases."));
+
+  return all_face_orientations;
 }
 
 
index cf562d5740801e0177961c83db51f3c710e69171..a434fe9b01907d10101f836527eb83f3ee806863 100644 (file)
@@ -1179,7 +1179,7 @@ namespace internal
         MappingInfoStorage<dim, dim, double, VectorizedDouble> temp_data;
         temp_data.copy_descriptor(my_data);
         FEEvaluationBaseData<dim, double, false, VectorizedDouble> eval(
-          {&shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0);
+          {&shape_info, nullptr, &temp_data, 0, 0});
 
         AlignedVector<VectorizedDouble> evaluation_data;
         eval.set_data_pointers(&evaluation_data, dim);
@@ -2091,9 +2091,9 @@ namespace internal
         MappingInfoStorage<dim - 1, dim, double, VectorizedDouble> temp_data;
         temp_data.copy_descriptor(my_data);
         FEEvaluationBaseData<dim, double, true, VectorizedDouble> eval_int(
-          {&shape_info, nullptr, &temp_data, 0, 0}, 0, true, 0);
+          {&shape_info, nullptr, &temp_data, 0, 0}, true);
         FEEvaluationBaseData<dim, double, true, VectorizedDouble> eval_ext(
-          {&shape_info, nullptr, &temp_data, 0, 0}, 0, false, 0);
+          {&shape_info, nullptr, &temp_data, 0, 0}, false);
 
         AlignedVector<VectorizedDouble> evaluation_data, evaluation_data_ext;
         eval_int.set_data_pointers(&evaluation_data, dim);
index 6706195df71f9f6b4f2d782623c3e2be075a45c0..4ce2186a4e2ae4883ffeb50d534dc12cef851e25 100644 (file)
@@ -120,14 +120,12 @@ test(const unsigned int n_refinements = 1)
                                                    VectorizedArrayType>::
                 template interpolate_quadrature<true, false>(
                   1,
+                  EvaluationFlags::values,
                   matrix_free.get_shape_info(),
                   phi.begin_values(),
                   temp.data(),
-                  false,
-                  false,
                   face);
 
-
               for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
                 {
                   const auto u_cell = temp[q];

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.