From: Maximilian Bergbauer <maximilian.bergbauer@tum.de>
Date: Thu, 29 Jun 2023 14:19:36 +0000 (+0200)
Subject: Fix ECL in read_write_operation_contiguous
X-Git-Tag: relicensing~767^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15527%2Fhead;p=dealii.git

Fix ECL in read_write_operation_contiguous
---

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index a171bcda11..cf06095a7d 100644
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -3370,6 +3370,11 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                                            *this->dof_info);
     }
 
+  const bool accesses_exterior_dofs =
+    this->dof_access_index ==
+      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+    this->is_interior_face() == false;
+
   // Case 2: contiguous indices which use reduced storage of indices and can
   // use vectorized load/store operations -> go to separate function
   if (this->cell != numbers::invalid_unsigned_int)
@@ -3380,9 +3385,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
       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))
+      if (accesses_exterior_dofs)
         {
           const std::array<unsigned int, VectorizedArrayType::size()> &cells =
             this->get_cell_ids();
@@ -3451,7 +3454,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                          internal::is_vectorizable<VectorType, Number>::value>
     vector_selector;
 
-  const bool use_vectorized_path = !(masking_is_active || has_hn_constraints);
+  const bool use_vectorized_path =
+    !(masking_is_active || has_hn_constraints || accesses_exterior_dofs);
 
   const std::size_t dofs_per_component = this->data->dofs_per_component_on_cell;
   std::array<VectorizedArrayType *, n_components> values_dofs;
@@ -3845,17 +3849,17 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
   Assert(this->cell != numbers::invalid_unsigned_int, ExcNotImplemented());
 
+  const bool accesses_exterior_dofs =
+    this->dof_access_index ==
+      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+    this->is_interior_face() == false;
+
   // Simple case: We have contiguous storage, so we can simply copy out the
   // data
-  if ((dof_info.index_storage_variants[ind][this->cell] ==
-         internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-           interleaved_contiguous &&
-       n_lanes == VectorizedArrayType::size()) &&
-      !(is_face &&
-        this->dof_access_index ==
-          internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-        this->is_interior_face() == false) &&
-      !(!is_face && !this->is_interior_face()))
+  if (dof_info.index_storage_variants[ind][this->cell] ==
+        internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+          interleaved_contiguous &&
+      n_lanes == VectorizedArrayType::size() && !accesses_exterior_dofs)
     {
       const unsigned int dof_index =
         dof_indices_cont[this->cell * VectorizedArrayType::size()] +
@@ -3887,12 +3891,6 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   const unsigned int n_filled_lanes =
     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) ||
-    (!is_face && !this->is_interior_face());
-
   if (vectors_sm[0] != nullptr)
     {
       const auto compute_vector_ptrs = [&](const unsigned int comp) {
@@ -3941,7 +3939,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
       };
 
       if (n_filled_lanes == VectorizedArrayType::size() &&
-          n_lanes == VectorizedArrayType::size() && !is_ecl)
+          n_lanes == VectorizedArrayType::size() && !accesses_exterior_dofs)
         {
           if (n_components == 1 || this->n_fe_components == 1)
             {
@@ -4018,7 +4016,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   // In the case with contiguous cell indices, we know that there are no
   // constraints and that the indices within each element are contiguous
   if (n_filled_lanes == VectorizedArrayType::size() &&
-      n_lanes == VectorizedArrayType::size() && !is_ecl)
+      n_lanes == VectorizedArrayType::size() && !accesses_exterior_dofs)
     {
       if (dof_info.index_storage_variants[ind][this->cell] ==
           internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
@@ -4107,56 +4105,109 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
       {
         for (unsigned int i = 0; i < dofs_per_component; ++i)
           operation.process_empty(values_dofs[comp][i]);
-        if (dof_info.index_storage_variants[ind][this->cell] ==
-            internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-              contiguous)
-          {
-            if (n_components == 1 || this->n_fe_components == 1)
-              {
-                for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                  if (mask[v] == true)
-                    for (unsigned int i = 0; i < dofs_per_component; ++i)
-                      operation.process_dof(dof_indices[v] + i,
-                                            *src[comp],
-                                            values_dofs[comp][i][v]);
-              }
-            else
-              {
-                for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                  if (mask[v] == true)
-                    for (unsigned int i = 0; i < dofs_per_component; ++i)
-                      operation.process_dof(dof_indices[v] + i +
-                                              comp * dofs_per_component,
-                                            *src[0],
-                                            values_dofs[comp][i][v]);
-              }
-          }
-        else
+        if (accesses_exterior_dofs)
           {
-            const unsigned int *offsets =
-              &dof_info.dof_indices_interleave_strides
-                 [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 || this->n_fe_components == 1)
-              for (unsigned int v = 0; v < n_filled_lanes; ++v)
+              if (mask[v] == true)
                 {
-                  if (mask[v] == true)
-                    for (unsigned int i = 0; i < dofs_per_component; ++i)
-                      operation.process_dof(dof_indices[v] + i * offsets[v],
-                                            *src[comp],
-                                            values_dofs[comp][i][v]);
+                  if (dof_info.index_storage_variants
+                        [ind][cells[v] / VectorizedArrayType::size()] ==
+                      internal::MatrixFreeFunctions::DoFInfo::
+                        IndexStorageVariants::contiguous)
+                    {
+                      if (n_components == 1 || this->n_fe_components == 1)
+                        {
+                          for (unsigned int i = 0; i < dofs_per_component; ++i)
+                            operation.process_dof(dof_indices[v] + i,
+                                                  *src[comp],
+                                                  values_dofs[comp][i][v]);
+                        }
+                      else
+                        {
+                          for (unsigned int i = 0; i < dofs_per_component; ++i)
+                            operation.process_dof(dof_indices[v] + i +
+                                                    comp * dofs_per_component,
+                                                  *src[0],
+                                                  values_dofs[comp][i][v]);
+                        }
+                    }
+                  else
+                    {
+                      const unsigned int offset =
+                        dof_info.dof_indices_interleave_strides[ind][cells[v]];
+                      AssertIndexRange(offset, VectorizedArrayType::size() + 1);
+                      if (n_components == 1 || this->n_fe_components == 1)
+                        {
+                          for (unsigned int i = 0; i < dofs_per_component; ++i)
+                            operation.process_dof(dof_indices[v] + i * offset,
+                                                  *src[comp],
+                                                  values_dofs[comp][i][v]);
+                        }
+                      else
+                        {
+                          for (unsigned int i = 0; i < dofs_per_component; ++i)
+                            operation.process_dof(
+                              dof_indices[v] +
+                                (i + comp * dofs_per_component) * offset,
+                              *src[0],
+                              values_dofs[comp][i][v]);
+                        }
+                    }
                 }
+          }
+        else
+          {
+            if (dof_info.index_storage_variants[ind][this->cell] ==
+                internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                  contiguous)
+              {
+                if (n_components == 1 || this->n_fe_components == 1)
+                  {
+                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                      if (mask[v] == true)
+                        for (unsigned int i = 0; i < dofs_per_component; ++i)
+                          operation.process_dof(dof_indices[v] + i,
+                                                *src[comp],
+                                                values_dofs[comp][i][v]);
+                  }
+                else
+                  {
+                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                      if (mask[v] == true)
+                        for (unsigned int i = 0; i < dofs_per_component; ++i)
+                          operation.process_dof(dof_indices[v] + i +
+                                                  comp * dofs_per_component,
+                                                *src[0],
+                                                values_dofs[comp][i][v]);
+                  }
+              }
             else
               {
+                const unsigned int *offsets =
+                  &dof_info.dof_indices_interleave_strides
+                     [ind][VectorizedArrayType::size() * this->cell];
                 for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                  if (mask[v] == true)
-                    for (unsigned int i = 0; i < dofs_per_component; ++i)
-                      operation.process_dof(dof_indices[v] +
-                                              (i + comp * dofs_per_component) *
-                                                offsets[v],
-                                            *src[0],
-                                            values_dofs[comp][i][v]);
+                  AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
+                if (n_components == 1 || this->n_fe_components == 1)
+                  for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                    {
+                      if (mask[v] == true)
+                        for (unsigned int i = 0; i < dofs_per_component; ++i)
+                          operation.process_dof(dof_indices[v] + i * offsets[v],
+                                                *src[comp],
+                                                values_dofs[comp][i][v]);
+                    }
+                else
+                  {
+                    for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                      if (mask[v] == true)
+                        for (unsigned int i = 0; i < dofs_per_component; ++i)
+                          operation.process_dof(
+                            dof_indices[v] +
+                              (i + comp * dofs_per_component) * offsets[v],
+                            *src[0],
+                            values_dofs[comp][i][v]);
+                  }
               }
           }
       }
diff --git a/tests/matrix_free/ecl_05.cc b/tests/matrix_free/ecl_05.cc
new file mode 100644
index 0000000000..a91d93c784
--- /dev/null
+++ b/tests/matrix_free/ecl_05.cc
@@ -0,0 +1,51 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "ecl_fe_q.h"
+
+
+// Compare the evaluation of a SIP Laplace operator with face- and element-
+// centric loops on a hypercube, hypershell, and hyperball for FE_Q.
+// Adapted from ecl_01.cc
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  mpi_initlog();
+  deallog.precision(10);
+
+  test<2, 2, 3, double, VectorizedArray<double>>(0);
+  test<2, 2, 3, double, VectorizedArray<double>>(1);
+  test<2, 2, 3, double, VectorizedArray<double>>(2);
+}
diff --git a/tests/matrix_free/ecl_05.with_p4est=true.mpirun=1.output b/tests/matrix_free/ecl_05.with_p4est=true.mpirun=1.output
new file mode 100644
index 0000000000..a1d107322b
--- /dev/null
+++ b/tests/matrix_free/ecl_05.with_p4est=true.mpirun=1.output
@@ -0,0 +1,31 @@
+
+DEAL::36.95643080
+Process #0
+Local range: [0, 25), global size: 25
+Vector data:
+7.0000000000e+00 -8.0000000000e+00 6.6666666667e+00 -1.3333333333e+00 1.4666666667e+01 -5.3333333333e+00 6.6666666667e-01 1.3333333333e+00 2.6666666667e+00 7.0000000000e+00 6.6666666667e+00 1.4666666667e+01 6.6666666667e-01 1.3333333333e+00 2.6666666667e+00 7.0000000000e+00 -8.0000000000e+00 1.4666666667e+01 -5.3333333333e+00 6.6666666667e-01 2.6666666667e+00 7.0000000000e+00 1.4666666667e+01 6.6666666667e-01 2.6666666667e+00 
+DEAL::36.95643080
+Process #0
+Local range: [0, 25), global size: 25
+Vector data:
+7.0000000000e+00 -8.0000000000e+00 6.6666666667e+00 -1.3333333333e+00 1.4666666667e+01 -5.3333333333e+00 6.6666666667e-01 1.3333333333e+00 2.6666666667e+00 7.0000000000e+00 6.6666666667e+00 1.4666666667e+01 6.6666666667e-01 1.3333333333e+00 2.6666666667e+00 7.0000000000e+00 -8.0000000000e+00 1.4666666667e+01 -5.3333333333e+00 6.6666666667e-01 2.6666666667e+00 7.0000000000e+00 1.4666666667e+01 6.6666666667e-01 2.6666666667e+00 
+DEAL::76.88408355
+Process #0
+Local range: [0, 200), global size: 200
+Vector data:
+9.8820350839e+00 7.5921669161e+00 -2.4753536069e-01 -2.7024722797e-01 2.1654342873e+00 1.5967559056e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 -4.0795733480e+00 -4.5435251233e+00 -1.1272821381e+00 -1.1870898951e+00 -9.3180326752e+00 -2.3474056733e+00 -3.4516833425e+00 -1.1198445186e+00 -8.5073996496e+00 -2.3065557654e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 -9.1039260325e-01 -4.0990803825e-01 -5.1939423684e+00 -1.6282583741e+00 2.4174012868e+00 6.5059888355e-01 1.9926091726e+00 2.3474832686e-01 1.0032897601e+01 -4.4055047995e-01 2.5256889344e+00 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 4.9401455209e+00 1.3282850173e+00 8.1338631783e+00 2.1420439988e+00 2.4174012868e+00 6.5059888355e-01 8.1338631783e+00 2.1420439988e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -9.1039260325e-01 -4.0990803825e-01 1.9926091726e+00 2.3474832686e-01 -3.4516833425e+00 -1.1198445186e+00 -5.1939423684e+00 -1.6282583741e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 9.8820350839e+00 -2.4753536069e-01 2.1654342873e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 -4.5435251233e+00 -1.1870898951e+00 -8.5073996496e+00 -2.3065557654e+00 -4.0795733480e+00 -1.1272821381e+00 -9.3180326752e+00 -2.3474056733e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 -4.5435251233e+00 -1.1870898951e+00 -9.3180326752e+00 -2.3474056733e+00 -3.4516833425e+00 -1.1198445186e+00 -8.5073996496e+00 -2.3065557654e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 -9.1039260325e-01 -4.0990803825e-01 -5.1939423684e+00 -1.6282583741e+00 2.4174012868e+00 6.5059888355e-01 1.9926091726e+00 2.3474832686e-01 1.0032897601e+01 -4.4055047995e-01 2.5256889344e+00 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 4.9401455209e+00 1.3282850173e+00 8.1338631783e+00 2.1420439988e+00 2.4174012868e+00 6.5059888355e-01 8.1338631783e+00 2.1420439988e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -9.1039260325e-01 -4.0990803825e-01 1.9926091726e+00 2.3474832686e-01 -3.4516833425e+00 -1.1198445186e+00 -5.1939423684e+00 -1.6282583741e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 -4.5435251233e+00 -1.1870898951e+00 -8.5073996496e+00 -2.3065557654e+00 -9.3180326752e+00 -2.3474056733e+00 
+DEAL::76.88408355
+Process #0
+Local range: [0, 200), global size: 200
+Vector data:
+9.8820350839e+00 7.5921669161e+00 -2.4753536069e-01 -2.7024722797e-01 2.1654342873e+00 1.5967559056e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 -4.0795733480e+00 -4.5435251233e+00 -1.1272821381e+00 -1.1870898951e+00 -9.3180326752e+00 -2.3474056733e+00 -3.4516833425e+00 -1.1198445186e+00 -8.5073996496e+00 -2.3065557654e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 -9.1039260325e-01 -4.0990803825e-01 -5.1939423684e+00 -1.6282583741e+00 2.4174012868e+00 6.5059888355e-01 1.9926091726e+00 2.3474832686e-01 1.0032897601e+01 -4.4055047995e-01 2.5256889344e+00 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 4.9401455209e+00 1.3282850173e+00 8.1338631783e+00 2.1420439988e+00 2.4174012868e+00 6.5059888355e-01 8.1338631783e+00 2.1420439988e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -9.1039260325e-01 -4.0990803825e-01 1.9926091726e+00 2.3474832686e-01 -3.4516833425e+00 -1.1198445186e+00 -5.1939423684e+00 -1.6282583741e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 9.8820350839e+00 -2.4753536069e-01 2.1654342873e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 -4.5435251233e+00 -1.1870898951e+00 -8.5073996496e+00 -2.3065557654e+00 -4.0795733480e+00 -1.1272821381e+00 -9.3180326752e+00 -2.3474056733e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 -4.5435251233e+00 -1.1870898951e+00 -9.3180326752e+00 -2.3474056733e+00 -3.4516833425e+00 -1.1198445186e+00 -8.5073996496e+00 -2.3065557654e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 -9.1039260325e-01 -4.0990803825e-01 -5.1939423684e+00 -1.6282583741e+00 2.4174012868e+00 6.5059888355e-01 1.9926091726e+00 2.3474832686e-01 1.0032897601e+01 -4.4055047995e-01 2.5256889344e+00 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 -3.5555293364e+00 6.4394814131e-02 -7.2756796289e-01 9.6123319432e+00 -4.6799564242e-01 2.4984783418e+00 4.9401455209e+00 1.3282850173e+00 8.1338631783e+00 2.1420439988e+00 2.4174012868e+00 6.5059888355e-01 8.1338631783e+00 2.1420439988e+00 -9.2223610976e+00 2.4809504401e-01 -2.3412456432e+00 -1.5845965268e+01 4.4377777377e-01 -3.7720426578e+00 3.2500096518e+00 -1.8049079881e-01 4.6688151672e-01 -7.1664026326e+00 7.3617621377e-02 -2.1853085730e+00 -9.1039260325e-01 -4.0990803825e-01 1.9926091726e+00 2.3474832686e-01 -3.4516833425e+00 -1.1198445186e+00 -5.1939423684e+00 -1.6282583741e+00 7.5921669161e+00 -2.7024722797e-01 1.5967559056e+00 1.1163432635e+01 -4.6965932229e-01 2.1388829043e+00 1.7158598359e+01 -5.0563990584e-01 3.8153850674e+00 -4.5435251233e+00 -1.1870898951e+00 -8.5073996496e+00 -2.3065557654e+00 -9.3180326752e+00 -2.3474056733e+00 
+DEAL::59.34388966
+Process #0
+Local range: [0, 89), global size: 89
+Vector data:
+-1.1986619572e+01 1.5004834009e+01 -1.4604443071e-01 -1.6361280608e-01 -3.2544086977e+00 3.9763571412e+00 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 -1.1986619572e+01 -1.4604443071e-01 -3.2544086977e+00 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 -1.4524926676e-01 4.0635119219e-01 -8.4386266216e-01 8.3824944782e-01 4.2547396644e-01 6.8707503578e-01 -1.4524926676e-01 -8.4386266216e-01 4.2547396644e-01 6.8707503578e-01 -6.5212979468e+00 -1.8020594332e-01 -1.3042595894e+01 -3.6041188664e-01 -2.0542935534e+00 -4.1085871068e+00 -8.4128040724e-02 -1.6825608145e-01 -5.5321803262e-01 -1.1064360652e+00 -1.1986619572e+01 -1.4604443071e-01 -1.3042595894e+01 -3.6041188664e-01 -3.2544086977e+00 -4.1085871068e+00 -1.4524926676e-01 -1.6825608145e-01 -8.4386266216e-01 -1.1064360652e+00 4.1958793807e-01 8.3917587615e-01 4.2458228995e-01 8.4916457991e-01 -8.4128040724e-02 -1.6825608145e-01 4.2458228995e-01 8.4916457991e-01 4.0635119219e-01 8.3917587615e-01 4.2547396644e-01 8.4916457991e-01 -1.4524926676e-01 -1.6825608145e-01 4.2547396644e-01 8.4916457991e-01 -6.5212979468e+00 -1.8020594332e-01 -2.0542935534e+00 -1.3042595894e+01 -3.6041188664e-01 -4.1085871068e+00 -1.1986619572e+01 -1.4604443071e-01 -3.2544086977e+00 -1.3042595894e+01 -3.6041188664e-01 -4.1085871068e+00 -5.5321803262e-01 -1.1064360652e+00 -8.4386266216e-01 -1.1064360652e+00 1.5004834009e+01 -1.6361280608e-01 -1.9357889961e+01 1.1562719635e+00 3.9763571412e+00 -3.8249201375e+00 8.3824944782e-01 6.8707503578e-01 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 6.8707503578e-01 
+DEAL::59.34388966
+Process #0
+Local range: [0, 89), global size: 89
+Vector data:
+-1.1986619572e+01 1.5004834009e+01 -1.4604443071e-01 -1.6361280608e-01 -3.2544086977e+00 3.9763571412e+00 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 -1.1986619572e+01 -1.4604443071e-01 -3.2544086977e+00 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 -1.4524926676e-01 4.0635119219e-01 -8.4386266216e-01 8.3824944782e-01 4.2547396644e-01 6.8707503578e-01 -1.4524926676e-01 -8.4386266216e-01 4.2547396644e-01 6.8707503578e-01 -6.5212979468e+00 -1.8020594332e-01 -1.3042595894e+01 -3.6041188664e-01 -2.0542935534e+00 -4.1085871068e+00 -8.4128040724e-02 -1.6825608145e-01 -5.5321803262e-01 -1.1064360652e+00 -1.1986619572e+01 -1.4604443071e-01 -1.3042595894e+01 -3.6041188664e-01 -3.2544086977e+00 -4.1085871068e+00 -1.4524926676e-01 -1.6825608145e-01 -8.4386266216e-01 -1.1064360652e+00 4.1958793807e-01 8.3917587615e-01 4.2458228995e-01 8.4916457991e-01 -8.4128040724e-02 -1.6825608145e-01 4.2458228995e-01 8.4916457991e-01 4.0635119219e-01 8.3917587615e-01 4.2547396644e-01 8.4916457991e-01 -1.4524926676e-01 -1.6825608145e-01 4.2547396644e-01 8.4916457991e-01 -6.5212979468e+00 -1.8020594332e-01 -2.0542935534e+00 -1.3042595894e+01 -3.6041188664e-01 -4.1085871068e+00 -1.1986619572e+01 -1.4604443071e-01 -3.2544086977e+00 -1.3042595894e+01 -3.6041188664e-01 -4.1085871068e+00 -5.5321803262e-01 -1.1064360652e+00 -8.4386266216e-01 -1.1064360652e+00 1.5004834009e+01 -1.6361280608e-01 -1.9357889961e+01 1.1562719635e+00 3.9763571412e+00 -3.8249201375e+00 8.3824944782e-01 6.8707503578e-01 -1.9357889961e+01 1.1562719635e+00 -3.8249201375e+00 6.8707503578e-01 
diff --git a/tests/matrix_free/ecl_fe_q.h b/tests/matrix_free/ecl_fe_q.h
new file mode 100644
index 0000000000..7a29312e68
--- /dev/null
+++ b/tests/matrix_free/ecl_fe_q.h
@@ -0,0 +1,346 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/tools.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+// Adapted from ecl.h
+
+template <int dim, typename Number>
+class CosineFunction : public Function<dim, Number>
+{
+public:
+  Number
+  value(const Point<dim> &p, const unsigned int component = 0) const
+  {
+    (void)component;
+    return std::cos(2 * numbers::PI * p[0]);
+  }
+};
+
+template <int dim,
+          int fe_degree,
+          int n_points                 = fe_degree + 1,
+          typename Number              = double,
+          typename VectorizedArrayType = VectorizedArray<Number>>
+void
+test(const unsigned int geometry      = 0,
+     const unsigned int n_refinements = 1,
+     const bool         print_vector  = true,
+     const MPI_Comm     comm          = MPI_COMM_SELF)
+{
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+
+  if (geometry == 0)
+    GridGenerator::hyper_cube(tria, 0.0, 1.0, true);
+  else if (geometry == 1)
+    GridGenerator::hyper_shell(tria, Point<dim>(), 0.5, 1.0);
+  else if (geometry == 2)
+    GridGenerator::hyper_ball(tria);
+  else
+    Assert(false, ExcNotImplemented());
+
+  tria.reset_all_manifolds();
+
+  if (false)
+    {
+      std::vector<dealii::GridTools::PeriodicFacePair<
+        typename dealii::Triangulation<dim>::cell_iterator>>
+        periodic_faces;
+
+      if (dim >= 1)
+        dealii::GridTools::collect_periodic_faces(
+          tria, 0, 1, 0, periodic_faces);
+
+      if (dim >= 2)
+        dealii::GridTools::collect_periodic_faces(
+          tria, 2, 3, 1, periodic_faces);
+
+      if (dim >= 3)
+        dealii::GridTools::collect_periodic_faces(
+          tria, 4, 5, 2, periodic_faces);
+
+      tria.add_periodicity(periodic_faces);
+    }
+
+  tria.refine_global(n_refinements);
+
+  FE_Q<dim>       fe(fe_degree);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  QGauss<1>     quad(n_points);
+
+  AffineConstraints<Number> constraint;
+
+  using MF = MatrixFree<dim, Number, VectorizedArrayType>;
+
+  typename MF::AdditionalData additional_data;
+  additional_data.mapping_update_flags = update_values | update_gradients;
+  additional_data.mapping_update_flags_inner_faces =
+    update_values | update_gradients;
+  additional_data.mapping_update_flags_boundary_faces =
+    update_values | update_gradients;
+  additional_data.tasks_parallel_scheme =
+    MF::AdditionalData::TasksParallelScheme::none;
+  additional_data.mapping_update_flags_faces_by_cells =
+    update_values | update_gradients;
+  additional_data.hold_all_faces_to_owned_cells = true;
+  additional_data.communicator_sm               = comm;
+
+  MF matrix_free;
+  matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data);
+
+  VectorType src, dst;
+
+  matrix_free.initialize_dof_vector(src);
+  matrix_free.initialize_dof_vector(dst);
+
+  FEEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType> phi(
+    matrix_free);
+  FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+    phi_m(matrix_free, true);
+  FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+    phi_p(matrix_free, false);
+
+
+  CosineFunction<dim, Number> function;
+  VectorTools::interpolate(dof_handler, function, src);
+
+  dst = 0.0;
+
+  /**
+   * Face-centric loop
+   */
+  matrix_free.template loop<VectorType, VectorType>(
+    [&](const auto &, auto &dst, const auto &src, const auto range) {
+      for (unsigned int cell = range.first; cell < range.second; ++cell)
+        {
+          phi.reinit(cell);
+          phi.read_dof_values(src);
+          phi.evaluate(EvaluationFlags::gradients);
+          for (unsigned int q = 0; q < phi.n_q_points; ++q)
+            phi.submit_gradient(phi.get_gradient(q), q);
+          phi.integrate(EvaluationFlags::gradients);
+          phi.distribute_local_to_global(dst);
+        }
+    },
+    [&](const auto &, auto &dst, const auto &src, const auto range) {
+      for (unsigned int face = range.first; face < range.second; ++face)
+        {
+          phi_m.reinit(face);
+          phi_p.reinit(face);
+
+          phi_m.read_dof_values(src);
+          phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+          phi_p.read_dof_values(src);
+          phi_p.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+          VectorizedArrayType sigmaF =
+            (std::abs(
+               (phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) +
+             std::abs(
+               (phi_m.normal_vector(0) * phi_p.inverse_jacobian(0))[dim - 1])) *
+            (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0));
+
+          for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+            {
+              VectorizedArrayType average_value =
+                (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5;
+              VectorizedArrayType average_valgrad =
+                phi_m.get_normal_derivative(q) + phi_p.get_normal_derivative(q);
+              average_valgrad =
+                average_value * 2. * sigmaF - average_valgrad * 0.5;
+              phi_m.submit_normal_derivative(-average_value, q);
+              phi_p.submit_normal_derivative(-average_value, q);
+              phi_m.submit_value(average_valgrad, q);
+              phi_p.submit_value(-average_valgrad, q);
+            }
+          phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+          phi_m.distribute_local_to_global(dst);
+          phi_p.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+          phi_p.distribute_local_to_global(dst);
+        }
+    },
+    [&](const auto &, auto &dst, const auto &src, const auto face_range) {
+      for (unsigned int face = face_range.first; face < face_range.second;
+           face++)
+        {
+          phi_m.reinit(face);
+          phi_m.read_dof_values(src);
+          phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+          VectorizedArrayType sigmaF =
+            std::abs(
+              (phi_m.normal_vector(0) * phi_m.inverse_jacobian(0))[dim - 1]) *
+            Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0;
+
+          for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+            {
+              VectorizedArrayType average_value = phi_m.get_value(q);
+              VectorizedArrayType average_valgrad =
+                -phi_m.get_normal_derivative(q);
+              average_valgrad += average_value * sigmaF * 2.0;
+              phi_m.submit_normal_derivative(-average_value, q);
+              phi_m.submit_value(average_valgrad, q);
+            }
+
+          phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+          phi_m.distribute_local_to_global(dst);
+        }
+    },
+    dst,
+    src,
+    false,
+    MF::DataAccessOnFaces::gradients,
+    MF::DataAccessOnFaces::gradients);
+
+  deallog << dst.l2_norm() << std::endl;
+
+  if (print_vector)
+    dst.print(deallog.get_file_stream(), 10);
+
+  dst = 0.0;
+
+  /**
+   * Element-centric loop
+   */
+  matrix_free.template loop_cell_centric<VectorType, VectorType>(
+    [&](const auto &, auto &dst, const auto &src, const auto range) {
+      for (unsigned int cell = range.first; cell < range.second; ++cell)
+        {
+          phi.reinit(cell);
+          phi.read_dof_values(src);
+          phi.evaluate(EvaluationFlags::gradients);
+          for (unsigned int q = 0; q < phi.n_q_points; ++q)
+            phi.submit_gradient(phi.get_gradient(q), q);
+          phi.integrate(EvaluationFlags::gradients);
+
+          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+               ++face)
+            {
+              const auto bids =
+                matrix_free.get_faces_by_cells_boundary_id(cell, face);
+
+              std::bitset<VectorizedArrayType::size()> mask;
+              for (unsigned int v = 0;
+                   v < matrix_free.n_active_entries_per_cell_batch(cell);
+                   ++v)
+                if (bids[v] == numbers::internal_face_boundary_id)
+                  mask[v] = true;
+
+              // boundary
+              {
+                phi_m.reinit(cell, face);
+                phi_m.read_dof_values(src);
+                phi_m.evaluate(EvaluationFlags::values |
+                               EvaluationFlags::gradients);
+                VectorizedArrayType sigmaF =
+                  std::abs((phi_m.normal_vector(0) *
+                            phi_m.inverse_jacobian(0))[dim - 1]) *
+                  Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0;
+
+                for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+                  {
+                    VectorizedArrayType average_value = phi_m.get_value(q);
+                    VectorizedArrayType average_valgrad =
+                      -phi_m.get_normal_derivative(q);
+                    average_valgrad += average_value * sigmaF * 2.0;
+                    phi_m.submit_normal_derivative(-average_value, q);
+                    phi_m.submit_value(average_valgrad, q);
+                  }
+
+                phi_m.integrate(EvaluationFlags::values |
+                                EvaluationFlags::gradients);
+
+                for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+                  for (unsigned int q = 0; q < phi.dofs_per_cell; ++q)
+                    phi.begin_dof_values()[q][v] +=
+                      mask[v] == false ? phi_m.begin_dof_values()[q][v] : 0.;
+              }
+
+              // interior face
+              {
+                phi_m.reinit(cell, face);
+                phi_p.reinit(cell, face);
+
+                phi_m.read_dof_values(src);
+                phi_m.evaluate(EvaluationFlags::values |
+                               EvaluationFlags::gradients);
+
+                phi_p.read_dof_values(src, 0, mask);
+                phi_p.evaluate(EvaluationFlags::values |
+                               EvaluationFlags::gradients);
+
+                VectorizedArrayType sigmaF =
+                  (std::abs((phi_m.normal_vector(0) *
+                             phi_m.inverse_jacobian(0))[dim - 1]) +
+                   std::abs((phi_m.normal_vector(0) *
+                             phi_p.inverse_jacobian(0))[dim - 1])) *
+                  (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0));
+
+                for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+                  {
+                    VectorizedArrayType average_value =
+                      (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5;
+                    VectorizedArrayType average_valgrad =
+                      phi_m.get_normal_derivative(q) +
+                      phi_p.get_normal_derivative(q);
+                    average_valgrad =
+                      average_value * 2. * sigmaF - average_valgrad * 0.5;
+                    phi_m.submit_normal_derivative(-average_value, q);
+                    phi_m.submit_value(average_valgrad, q);
+                  }
+                phi_m.integrate(EvaluationFlags::values |
+                                EvaluationFlags::gradients);
+
+                for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+                  for (unsigned int q = 0; q < phi.dofs_per_cell; ++q)
+                    phi.begin_dof_values()[q][v] +=
+                      mask[v] == true ? phi_m.begin_dof_values()[q][v] : 0.;
+              }
+            }
+
+          phi.distribute_local_to_global(dst);
+        }
+    },
+    dst,
+    src,
+    false,
+    MF::DataAccessOnFaces::gradients);
+
+  deallog << dst.l2_norm() << std::endl;
+
+  if (print_vector)
+    dst.print(deallog.get_file_stream(), 10);
+}