Introduce is_cell_neighbor flag to select correct code path

include/deal.II/matrix_free/fe_evaluation.h

@@ -444,7 +444,9 @@ class FEEvaluationBaseData
 
   /**
    * Flag holding information whether a face is an interior or exterior face
-   * according to the defined direction of the normal.  Not used for cells.
+   * according to the defined direction of the normal. For cells it defines if
+   * the dof values should be read from the actual cell corresponding to the
+   * interior face or the neighboring cell corresponding to the exterior face
    */
   bool is_interior_face;
 
@@ -3964,10 +3966,13 @@ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
   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)
+  const bool is_cell_neighbor = !is_face && !this->is_interior_face;
+
+  if ((is_face &&
+       this->dof_access_index ==
+         internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+       this->is_interior_face == false) ||
+      is_cell_neighbor)
     {
       // cell-based face-loop: plus face
       for (unsigned int i = 0; i < v_len; ++i)
@@ -4920,6 +4925,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   const std::vector<unsigned int> &dof_indices_cont =
     this->dof_info->dof_indices_contiguous[ind];
 
+  const bool is_cell_neighbor = !is_face && !this->is_interior_face;
+
   // Simple case: We have contiguous storage, so we can simply copy out the
   // data
   if ((this->dof_info->index_storage_variants[ind][this->cell] ==
@@ -4929,7 +4936,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
       !(is_face &&
         this->dof_access_index ==
           internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-        this->is_interior_face == false))
+        this->is_interior_face == false) &&
+      (is_cell_neighbor == false))
     {
       const unsigned int dof_index =
         dof_indices_cont[this->cell * VectorizedArrayType::size()] +
@@ -4963,9 +4971,10 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
 
   const bool is_ecl =
-    this->dof_access_index ==
-      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-    this->is_interior_face == false;
+    (this->dof_access_index ==
+       internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+     this->is_interior_face == false) ||
+    is_cell_neighbor;
 
   if (vectors_sm[0] != nullptr)
     {