Select general path for exterior cells with constraints

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 077067d7a697692e667cdcd7e85f6c0d63413c28..74257d040a23dea10333f5ad7ecc2a641b882949 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -3377,11 +3377,40 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
       AssertIndexRange(
         this->cell,
         dof_info.index_storage_variants[this->dof_access_index].size());
-      if (dof_info.index_storage_variants
-            [is_face ? this->dof_access_index :
-                       internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-            [this->cell] >= internal::MatrixFreeFunctions::DoFInfo::
-                              IndexStorageVariants::contiguous)
+
+      bool is_contiguous = true;
+      // check if exterior cells are not contiguous (ECL case)
+      if (is_face && !this->interior_face &&
+          (this->dof_access_index ==
+           internal::MatrixFreeFunctions::DoFInfo::dof_access_cell))
+        {
+          const std::array<unsigned int, VectorizedArrayType::size()> &cells =
+            this->get_cell_ids();
+          const unsigned int n_filled_lanes =
+            dof_info.n_vectorization_lanes_filled
+              [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+              [this->cell];
+          // we have to check all filled lanes which are active in the mask
+          for (unsigned int v = 0; v < n_filled_lanes; ++v)
+            if (mask[v] == true &&
+                dof_info.index_storage_variants
+                    [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+                    [cells[v] / VectorizedArrayType::size()] <
+                  internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                    contiguous)
+              is_contiguous = false;
+        } // or if cell/face batch is not contiguous
+      else if (dof_info.index_storage_variants
+                 [is_face ?
+                    this->dof_access_index :
+                    internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+                 [this->cell] < internal::MatrixFreeFunctions::DoFInfo::
+                                  IndexStorageVariants::contiguous)
+        {
+          is_contiguous = false;
+        }
+
+      if (is_contiguous)
         {
           read_write_operation_contiguous(operation, src, src_sm, mask);
           return;