]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Generalize MatrixFreeTools::compute_diagonal() and MatrixFreeTools::compute_matrix() 17029/head
authorPeter Munch <peterrmuench@gmail.com>
Wed, 8 May 2024 22:26:04 +0000 (00:26 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Fri, 7 Jun 2024 12:25:09 +0000 (14:25 +0200)
doc/news/changes/minor/20240516Munch [new file with mode: 0644]
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/tools.h
tests/matrix_free/compute_diagonal_08.cc [new file with mode: 0644]
tests/matrix_free/compute_diagonal_08.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20240516Munch b/doc/news/changes/minor/20240516Munch
new file mode 100644 (file)
index 0000000..f141495
--- /dev/null
@@ -0,0 +1,4 @@
+New: MatrixFreeTools::compute_diagonal() and MatrixFreeTools::compute_matrix()
+can now handle face and boundary integrals needed for DG.
+<br>
+(Peter Munch, Magdalena Schreter-Fleischhacker, 2024/05/16)
index 6dae7061118b7c11bf8f06a5a553e326fca8c664..adfce0cb6dfa26b829b0138ae653a20f1770a589 100644 (file)
@@ -2744,6 +2744,20 @@ namespace internal
       return vec[component];
     }
   };
+
+  template <typename VectorType, std::size_t N>
+  struct BlockVectorSelector<std::array<VectorType *, N>, false>
+  {
+    using BaseVectorType = VectorType;
+
+    static BaseVectorType *
+    get_vector_component(std::array<VectorType *, N> &vec,
+                         const unsigned int           component)
+    {
+      AssertIndexRange(component, vec.size());
+      return vec[component];
+    }
+  };
 } // namespace internal
 
 
index 73f77159f694603d4b581827b71fabc336101427..92e3f83166ea88537f3292653e6c56f80f6a83e0 100644 (file)
@@ -32,6 +32,31 @@ DEAL_II_NAMESPACE_OPEN
  */
 namespace MatrixFreeTools
 {
+  namespace internal
+  {
+    template <int dim, typename Number, bool is_face_>
+    class ComputeMatrixScratchData
+    {
+    public:
+      using FEEvalType = FEEvaluationData<dim, Number, is_face_>;
+
+      std::vector<unsigned int> dof_numbers;
+      std::vector<unsigned int> quad_numbers;
+      std::vector<unsigned int> first_selected_components;
+      std::vector<unsigned int> batch_type;
+      static const bool         is_face = is_face_;
+
+      std::function<std::vector<std::unique_ptr<FEEvalType>>(
+        const std::pair<unsigned int, unsigned int> &)>
+        op_create;
+      std::function<void(std::vector<std::unique_ptr<FEEvalType>> &,
+                         const unsigned int)>
+        op_reinit;
+      std::function<void(std::vector<std::unique_ptr<FEEvalType>> &)>
+        op_compute;
+    };
+  } // namespace internal
+
   /**
    * Modify @p additional_data so that cells are categorized
    * according to their boundary IDs, making face integrals in the case of
@@ -42,9 +67,11 @@ namespace MatrixFreeTools
   categorize_by_boundary_ids(const Triangulation<dim> &tria,
                              AdditionalData           &additional_data);
 
+
+
   /**
    * Compute the diagonal of a linear operator (@p diagonal_global), given
-   * @p matrix_free and the local cell integral operation @p local_vmult. The
+   * @p matrix_free and the local cell integral operation @p cell_operation. The
    * vector is initialized to the right size in the function.
    *
    * The parameters @p dof_no, @p quad_no, and @p first_selected_component are
@@ -66,11 +93,14 @@ namespace MatrixFreeTools
                                           n_q_points_1d,
                                           n_components,
                                           Number,
-                                          VectorizedArrayType> &)> &local_vmult,
-    const unsigned int                                              dof_no  = 0,
-    const unsigned int                                              quad_no = 0,
+                                          VectorizedArrayType> &)>
+                      &cell_operation,
+    const unsigned int dof_no                   = 0,
+    const unsigned int quad_no                  = 0,
     const unsigned int first_selected_component = 0);
 
+
+
   /**
    * Same as above but with a class and a function pointer.
    */
@@ -98,9 +128,112 @@ namespace MatrixFreeTools
     const unsigned int first_selected_component = 0);
 
 
+
+  /**
+   * Compute the diagonal of a linear operator (@p diagonal_global), given
+   * @p matrix_free and the local cell integral operation @p cell_operation,
+   * interior face integral operation @p face_operation, and boundary face
+   * integral operation @p boundary_operation. The
+   * vector is initialized to the right size in the function.
+   *
+   * The parameters @p dof_no, @p quad_no, and @p first_selected_component are
+   * passed to the constructor of the FEEvaluation that is internally set up.
+   */
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename VectorType>
+  void
+  compute_diagonal(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    VectorType                                         &diagonal_global,
+    const std::function<void(FEEvaluation<dim,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          n_components,
+                                          Number,
+                                          VectorizedArrayType> &)>
+      &cell_operation,
+    const std::function<void(FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &,
+                             FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &)>
+      &face_operation,
+    const std::function<void(FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &)>
+                      &boundary_operation,
+    const unsigned int dof_no                   = 0,
+    const unsigned int quad_no                  = 0,
+    const unsigned int first_selected_component = 0);
+
+
+
+  /**
+   * Same as above but with a class and a function pointer.
+   */
+  template <typename CLASS,
+            int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename VectorType>
+  void
+  compute_diagonal(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    VectorType                                         &diagonal_global,
+    void (CLASS::*cell_operation)(FEEvaluation<dim,
+                                               fe_degree,
+                                               n_q_points_1d,
+                                               n_components,
+                                               Number,
+                                               VectorizedArrayType> &) const,
+    void (CLASS::*face_operation)(FEFaceEvaluation<dim,
+                                                   fe_degree,
+                                                   n_q_points_1d,
+                                                   n_components,
+                                                   Number,
+                                                   VectorizedArrayType> &,
+                                  FEFaceEvaluation<dim,
+                                                   fe_degree,
+                                                   n_q_points_1d,
+                                                   n_components,
+                                                   Number,
+                                                   VectorizedArrayType> &)
+      const,
+    void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
+                                                       fe_degree,
+                                                       n_q_points_1d,
+                                                       n_components,
+                                                       Number,
+                                                       VectorizedArrayType> &)
+      const,
+    const CLASS       *owning_class,
+    const unsigned int dof_no                   = 0,
+    const unsigned int quad_no                  = 0,
+    const unsigned int first_selected_component = 0);
+
+
+
   /**
    * Compute the matrix representation of a linear operator (@p matrix), given
-   * @p matrix_free and the local cell integral operation @p local_vmult.
+   * @p matrix_free and the local cell integral operation @p cell_operation.
    * Constrained entries on the diagonal are set to one.
    *
    * The parameters @p dof_no, @p quad_no, and @p first_selected_component are
@@ -115,19 +248,22 @@ namespace MatrixFreeTools
             typename MatrixType>
   void
   compute_matrix(
-    const MatrixFree<dim, Number, VectorizedArrayType>             &matrix_free,
-    const AffineConstraints<Number>                                &constraints,
-    MatrixType                                                     &matrix,
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const AffineConstraints<Number>                    &constraints,
+    MatrixType                                         &matrix,
     const std::function<void(FEEvaluation<dim,
                                           fe_degree,
                                           n_q_points_1d,
                                           n_components,
                                           Number,
-                                          VectorizedArrayType> &)> &local_vmult,
-    const unsigned int                                              dof_no  = 0,
-    const unsigned int                                              quad_no = 0,
+                                          VectorizedArrayType> &)>
+                      &cell_operation,
+    const unsigned int dof_no                   = 0,
+    const unsigned int quad_no                  = 0,
     const unsigned int first_selected_component = 0);
 
+
+
   /**
    * Same as above but with a class and a function pointer.
    */
@@ -156,6 +292,133 @@ namespace MatrixFreeTools
     const unsigned int first_selected_component = 0);
 
 
+  namespace internal
+  {
+    /**
+     * Compute the matrix representation of a linear operator (@p matrix), given
+     * @p matrix_free and the local cell, face and boundary integral operation.
+     */
+    template <int dim,
+              typename Number,
+              typename VectorizedArrayType,
+              typename MatrixType>
+    void
+    compute_matrix(
+      const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+      const AffineConstraints<Number>                    &constraints,
+      const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
+        &cell_operation,
+      const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+        &face_operation,
+      const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+                 &boundary_operation,
+      MatrixType &matrix);
+  } // namespace internal
+
+
+
+  /**
+   * Compute the matrix representation of a linear operator (@p matrix), given
+   * @p matrix_free and the local cell integral operation @p cell_operation,
+   * interior face integral operation @p face_operation, and boundary face
+   * integal operation @p boundary_operation.
+   *
+   * The parameters @p dof_no, @p quad_no, and @p first_selected_component are
+   * passed to the constructor of the FEEvaluation that is internally set up.
+   */
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename MatrixType>
+  void
+  compute_matrix(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const AffineConstraints<Number>                    &constraints,
+    MatrixType                                         &matrix,
+    const std::function<void(FEEvaluation<dim,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          n_components,
+                                          Number,
+                                          VectorizedArrayType> &)>
+      &cell_operation,
+    const std::function<void(FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &,
+                             FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &)>
+      &face_operation,
+    const std::function<void(FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &)>
+                      &boundary_operation,
+    const unsigned int dof_no                   = 0,
+    const unsigned int quad_no                  = 0,
+    const unsigned int first_selected_component = 0);
+
+
+
+  /**
+   * Same as above but with a class and a function pointer.
+   */
+  template <typename CLASS,
+            int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename MatrixType>
+  void
+  compute_matrix(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const AffineConstraints<Number>                    &constraints,
+    MatrixType                                         &matrix,
+    void (CLASS::*cell_operation)(FEEvaluation<dim,
+                                               fe_degree,
+                                               n_q_points_1d,
+                                               n_components,
+                                               Number,
+                                               VectorizedArrayType> &) const,
+    void (CLASS::*face_operation)(FEFaceEvaluation<dim,
+                                                   fe_degree,
+                                                   n_q_points_1d,
+                                                   n_components,
+                                                   Number,
+                                                   VectorizedArrayType> &,
+                                  FEFaceEvaluation<dim,
+                                                   fe_degree,
+                                                   n_q_points_1d,
+                                                   n_components,
+                                                   Number,
+                                                   VectorizedArrayType> &)
+      const,
+    void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
+                                                       fe_degree,
+                                                       n_q_points_1d,
+                                                       n_components,
+                                                       Number,
+                                                       VectorizedArrayType> &)
+      const,
+    const CLASS       *owning_class,
+    const unsigned int dof_no                   = 0,
+    const unsigned int quad_no                  = 0,
+    const unsigned int first_selected_component = 0);
+
+
 
   /**
    * A wrapper around MatrixFree to help users to deal with DoFHandler
@@ -430,16 +693,16 @@ namespace MatrixFreeTools
       std::vector<std::pair<unsigned int, unsigned int>> inverse_lookup_origins;
     };
 
-    template <int dim,
-              int fe_degree,
-              int n_q_points_1d,
-              int n_components,
-              typename Number,
-              typename VectorizedArrayType>
+    template <typename FEEvaluationType, bool is_face>
     class ComputeDiagonalHelper
     {
     public:
-      static const unsigned int n_lanes = VectorizedArrayType::size();
+      using Number              = typename FEEvaluationType::number_type;
+      using VectorizedArrayType = typename FEEvaluationType::NumberType;
+
+      static const unsigned int dim          = FEEvaluationType::dimension;
+      static const unsigned int n_components = FEEvaluationType::n_components;
+      static const unsigned int n_lanes      = VectorizedArrayType::size();
 
       ComputeDiagonalHelper()
         : phi(nullptr)
@@ -452,12 +715,7 @@ namespace MatrixFreeTools
       {}
 
       void
-      initialize(FEEvaluation<dim,
-                              fe_degree,
-                              n_q_points_1d,
-                              n_components,
-                              Number,
-                              VectorizedArrayType> &phi)
+      initialize(FEEvaluationType &phi)
       {
         // if we are in hp mode and the number of unknowns changed, we must
         // clear the map of entries
@@ -486,8 +744,9 @@ namespace MatrixFreeTools
         const unsigned int first_selected_component =
           n_fe_components == 1 ? 0 : phi->get_first_selected_component();
 
-        const unsigned int n_lanes_filled =
-          matrix_free.n_active_entries_per_cell_batch(cell);
+        this->n_lanes_filled =
+          is_face ? matrix_free.n_active_entries_per_face_batch(cell) :
+                    matrix_free.n_active_entries_per_cell_batch(cell);
 
         // STEP 2: setup CSR storage of transposed locally-relevant
         //   constraint matrix
@@ -500,13 +759,19 @@ namespace MatrixFreeTools
         std::vector<unsigned int>  constraint_position;
         std::vector<unsigned char> is_constrained_hn;
 
+        const std::array<unsigned int, n_lanes> &cells =
+          this->phi->get_cell_ids();
+
         for (unsigned int v = 0; v < n_lanes_filled; ++v)
           {
+            Assert(cells[v] != numbers::invalid_unsigned_int,
+                   ExcInternalError());
+
             const unsigned int *dof_indices;
             unsigned int        index_indicators, next_index_indicators;
 
             const unsigned int start =
-              (cell * n_lanes + v) * n_fe_components + first_selected_component;
+              cells[v] * n_fe_components + first_selected_component;
             dof_indices =
               dof_info.dof_indices.data() + dof_info.row_starts[start].first;
             index_indicators      = dof_info.row_starts[start].second;
@@ -625,7 +890,7 @@ namespace MatrixFreeTools
                   [phi->get_active_fe_index()][first_selected_component])
               {
                 const auto mask =
-                  dof_info.hanging_node_constraint_masks[cell * n_lanes + v];
+                  dof_info.hanging_node_constraint_masks[cells[v]];
 
                 // cell has hanging nodes
                 if (mask != dealii::internal::MatrixFreeFunctions::
@@ -787,6 +1052,42 @@ namespace MatrixFreeTools
             c_pools[v].row_lid_to_gid.size() *
               (n_fe_components == 1 ? n_components : 1),
             Number(0.0));
+
+        // check if fast path can be taken via FEEvaluation
+        bool use_fast_path = true;
+
+        for (unsigned int v = 0; v < n_lanes_filled; ++v)
+          {
+            auto &c_pool = c_pools[v];
+
+            for (unsigned int i = 0; i < c_pool.row.size() - 1; ++i)
+              {
+                if ((c_pool.row[i + 1] - c_pool.row[i]) > 1)
+                  {
+                    use_fast_path = false;
+                    break;
+                  }
+                else if (((c_pool.row[i + 1] - c_pool.row[i]) == 1) &&
+                         (c_pool.val[c_pool.row[i]] != 1.0))
+                  {
+                    use_fast_path = false;
+                    break;
+                  }
+              }
+
+            if (use_fast_path == false)
+              break;
+          }
+
+        if (use_fast_path)
+          {
+            temp_values.resize(phi->dofs_per_cell);
+            return;
+          }
+        else
+          {
+            temp_values.clear();
+          }
       }
 
       void
@@ -866,9 +1167,23 @@ namespace MatrixFreeTools
         dof_values[i] = Number(1);
       }
 
+      void
+      zero_basis_vector()
+      {
+        VectorizedArrayType *dof_values = phi->begin_dof_values();
+        for (const unsigned int j : phi->dof_indices())
+          dof_values[j] = VectorizedArrayType();
+      }
+
       void
       submit()
       {
+        if (!temp_values.empty())
+          {
+            temp_values[i] = phi->begin_dof_values()[i];
+            return;
+          }
+
         // if we have a block vector with components with the same DoFHandler,
         // we need to figure out which component and which DoF within the
         // component are we currently considering
@@ -881,10 +1196,7 @@ namespace MatrixFreeTools
 
         // apply local constraint matrix from left and from right:
         // loop over all rows of transposed constrained matrix
-        for (unsigned int v = 0;
-             v < phi->get_matrix_free().n_active_entries_per_cell_batch(
-                   phi->get_current_cell_index());
-             ++v)
+        for (unsigned int v = 0; v < n_lanes_filled; ++v)
           {
             const auto &c_pool = c_pools[v];
 
@@ -913,14 +1225,21 @@ namespace MatrixFreeTools
       distribute_local_to_global(
         std::array<VectorType *, n_components> &diagonal_global)
       {
+        if (!temp_values.empty())
+          {
+            for (unsigned int j = 0; j < temp_values.size(); ++j)
+              phi->begin_dof_values()[j] = temp_values[j];
+
+            phi->distribute_local_to_global(diagonal_global);
+
+            return;
+          }
+
         // STEP 4: assembly results: add into global vector
         const unsigned int n_fe_components =
           phi->get_dof_info().start_components.back();
 
-        for (unsigned int v = 0;
-             v < phi->get_matrix_free().n_active_entries_per_cell_batch(
-                   phi->get_current_cell_index());
-             ++v)
+        for (unsigned int v = 0; v < n_lanes_filled; ++v)
           // if we have a block vector with components with the same
           // DoFHandler, we need to loop over all components manually and
           // need to apply the correct shift
@@ -937,18 +1256,21 @@ namespace MatrixFreeTools
                   [v][j + comp * c_pools[v].row_lid_to_gid.size()]);
       }
 
+      bool
+      use_fast_path() const
+      {
+        return !temp_values.empty();
+      }
+
     private:
-      FEEvaluation<dim,
-                   fe_degree,
-                   n_q_points_1d,
-                   n_components,
-                   Number,
-                   VectorizedArrayType> *phi;
+      FEEvaluationType *phi;
 
       unsigned int dofs_per_component;
 
       unsigned int i;
 
+      unsigned int n_lanes_filled;
+
       std::array<internal::LocalCSR<Number>, n_lanes> c_pools;
 
       // local storage: buffer so that we access the global vector once
@@ -962,9 +1284,52 @@ namespace MatrixFreeTools
         locally_relevant_constraints_hn_map;
 
       // scratch array
+      AlignedVector<VectorizedArrayType> temp_values;
       AlignedVector<VectorizedArrayType> values_dofs;
     };
 
+    template <bool is_face,
+              int  dim,
+              typename Number,
+              typename VectorizedArrayType>
+    bool
+    is_fe_nothing(
+      const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+      const std::pair<unsigned int, unsigned int>        &range,
+      const unsigned int                                  dof_no,
+      const unsigned int                                  quad_no,
+      const unsigned int first_selected_component,
+      const unsigned int fe_degree,
+      const unsigned int n_q_points_1d,
+      const bool         is_interior_face = true)
+    {
+      const unsigned int static_n_q_points =
+        is_face ? Utilities::pow(n_q_points_1d, dim - 1) :
+                  Utilities::pow(n_q_points_1d, dim);
+
+      unsigned int active_fe_index = 0;
+      if (!is_face)
+        active_fe_index = matrix_free.get_cell_active_fe_index(range);
+      else if (is_interior_face)
+        active_fe_index = matrix_free.get_face_range_category(range).first;
+      else
+        active_fe_index = matrix_free.get_face_range_category(range).second;
+
+      const auto init_data = dealii::internal::
+        extract_initialization_data<is_face, dim, Number, VectorizedArrayType>(
+          matrix_free,
+          dof_no,
+          first_selected_component,
+          quad_no,
+          fe_degree,
+          static_n_q_points,
+          active_fe_index,
+          numbers::invalid_unsigned_int /*active_quad_index*/,
+          numbers::invalid_unsigned_int /*face_type*/);
+
+      return init_data.shape_info->dofs_per_component_on_cell == 0;
+    }
+
   } // namespace internal
 
   template <int dim,
@@ -983,9 +1348,106 @@ namespace MatrixFreeTools
                                           n_q_points_1d,
                                           n_components,
                                           Number,
-                                          VectorizedArrayType> &)> &local_vmult,
-    const unsigned int                                              dof_no,
-    const unsigned int                                              quad_no,
+                                          VectorizedArrayType> &)>
+                      &cell_operation,
+    const unsigned int dof_no,
+    const unsigned int quad_no,
+    const unsigned int first_selected_component)
+  {
+    compute_diagonal<dim,
+                     fe_degree,
+                     n_q_points_1d,
+                     n_components,
+                     Number,
+                     VectorizedArrayType,
+                     VectorType>(matrix_free,
+                                 diagonal_global,
+                                 cell_operation,
+                                 {},
+                                 {},
+                                 dof_no,
+                                 quad_no,
+                                 first_selected_component);
+  }
+
+  template <typename CLASS,
+            int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename VectorType>
+  void
+  compute_diagonal(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    VectorType                                         &diagonal_global,
+    void (CLASS::*cell_operation)(FEEvaluation<dim,
+                                               fe_degree,
+                                               n_q_points_1d,
+                                               n_components,
+                                               Number,
+                                               VectorizedArrayType> &) const,
+    const CLASS       *owning_class,
+    const unsigned int dof_no,
+    const unsigned int quad_no,
+    const unsigned int first_selected_component)
+  {
+    compute_diagonal<dim,
+                     fe_degree,
+                     n_q_points_1d,
+                     n_components,
+                     Number,
+                     VectorizedArrayType,
+                     VectorType>(
+      matrix_free,
+      diagonal_global,
+      [&](auto &phi) { (owning_class->*cell_operation)(phi); },
+      dof_no,
+      quad_no,
+      first_selected_component);
+  }
+
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename VectorType>
+  void
+  compute_diagonal(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    VectorType                                         &diagonal_global,
+    const std::function<void(FEEvaluation<dim,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          n_components,
+                                          Number,
+                                          VectorizedArrayType> &)>
+      &cell_operation,
+    const std::function<void(FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &,
+                             FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &)>
+      &face_operation,
+    const std::function<void(FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &)>
+                      &boundary_operation,
+    const unsigned int dof_no,
+    const unsigned int quad_no,
     const unsigned int first_selected_component)
   {
     int dummy = 0;
@@ -1023,19 +1485,46 @@ namespace MatrixFreeTools
           *diagonal_global_components[0], matrix_free, dof_info);
       }
 
-    using Helper = internal::ComputeDiagonalHelper<dim,
+    using Helper =
+      internal::ComputeDiagonalHelper<FEEvaluation<dim,
                                                    fe_degree,
                                                    n_q_points_1d,
                                                    n_components,
                                                    Number,
-                                                   VectorizedArrayType>;
-
-    Threads::ThreadLocalStorage<Helper> scratch_data;
-    matrix_free.template cell_loop<VectorType, int>(
+                                                   VectorizedArrayType>,
+                                      false>;
+
+    using HelperFace =
+      internal::ComputeDiagonalHelper<FEFaceEvaluation<dim,
+                                                       fe_degree,
+                                                       n_q_points_1d,
+                                                       n_components,
+                                                       Number,
+                                                       VectorizedArrayType>,
+                                      true>;
+
+    Threads::ThreadLocalStorage<Helper>     scratch_data;
+    Threads::ThreadLocalStorage<HelperFace> scratch_data_m;
+    Threads::ThreadLocalStorage<HelperFace> scratch_data_p;
+
+    const auto cell_operation_wrapped =
       [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
           VectorType &,
           const int &,
-          const std::pair<unsigned int, unsigned int> &range) mutable {
+          const std::pair<unsigned int, unsigned int> &range) {
+        if (!cell_operation)
+          return;
+
+        // shortcut for FE_Nothing cells
+        if (internal::is_fe_nothing<false>(matrix_free,
+                                           range,
+                                           dof_no,
+                                           quad_no,
+                                           first_selected_component,
+                                           fe_degree,
+                                           n_q_points_1d))
+          return;
+
         Helper &helper = scratch_data.get();
 
         FEEvaluation<dim,
@@ -1054,16 +1543,165 @@ namespace MatrixFreeTools
             for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
               {
                 helper.prepare_basis_vector(i);
-                local_vmult(phi);
+                cell_operation(phi);
                 helper.submit();
               }
 
             helper.distribute_local_to_global(diagonal_global_components);
           }
-      },
-      diagonal_global,
-      dummy,
-      false);
+      };
+
+    const auto face_operation_wrapped =
+      [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+          VectorType &,
+          const int &,
+          const std::pair<unsigned int, unsigned int> &range) {
+        if (!face_operation)
+          return;
+
+        // shortcut for faces containing purely FE_Nothing
+        if (internal::is_fe_nothing<true>(matrix_free,
+                                          range,
+                                          dof_no,
+                                          quad_no,
+                                          first_selected_component,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          true) ||
+            internal::is_fe_nothing<true>(matrix_free,
+                                          range,
+                                          dof_no,
+                                          quad_no,
+                                          first_selected_component,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          false))
+          return;
+
+        HelperFace &helper_m = scratch_data_m.get();
+        HelperFace &helper_p = scratch_data_p.get();
+
+        FEFaceEvaluation<dim,
+                         fe_degree,
+                         n_q_points_1d,
+                         n_components,
+                         Number,
+                         VectorizedArrayType>
+          phi_m(matrix_free,
+                range,
+                true,
+                dof_no,
+                quad_no,
+                first_selected_component);
+        FEFaceEvaluation<dim,
+                         fe_degree,
+                         n_q_points_1d,
+                         n_components,
+                         Number,
+                         VectorizedArrayType>
+          phi_p(matrix_free,
+                range,
+                false,
+                dof_no,
+                quad_no,
+                first_selected_component);
+
+        helper_m.initialize(phi_m);
+        helper_p.initialize(phi_p);
+
+        for (unsigned int face = range.first; face < range.second; ++face)
+          {
+            helper_m.reinit(face);
+            helper_p.reinit(face);
+
+            // make check only if both adjacent cells have DoFs
+            Assert(helper_m.use_fast_path() && helper_p.use_fast_path(),
+                   ExcNotImplemented());
+
+            for (unsigned int i = 0; i < phi_m.dofs_per_cell; ++i)
+              {
+                helper_m.prepare_basis_vector(i);
+                helper_p.zero_basis_vector();
+                face_operation(phi_m, phi_p);
+                helper_m.submit();
+              }
+
+            helper_m.distribute_local_to_global(diagonal_global_components);
+
+            for (unsigned int i = 0; i < phi_p.dofs_per_cell; ++i)
+              {
+                helper_m.zero_basis_vector();
+                helper_p.prepare_basis_vector(i);
+                face_operation(phi_m, phi_p);
+                helper_p.submit();
+              }
+
+            helper_p.distribute_local_to_global(diagonal_global_components);
+          }
+      };
+
+    const auto boundary_operation_wrapped =
+      [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+          VectorType &,
+          const int &,
+          const std::pair<unsigned int, unsigned int> &range) {
+        if (!boundary_operation)
+          return;
+
+        // shortcut for faces containing purely FE_Nothing
+        if (internal::is_fe_nothing<true>(matrix_free,
+                                          range,
+                                          dof_no,
+                                          quad_no,
+                                          first_selected_component,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          true))
+          return;
+
+        HelperFace &helper = scratch_data_m.get();
+
+        FEFaceEvaluation<dim,
+                         fe_degree,
+                         n_q_points_1d,
+                         n_components,
+                         Number,
+                         VectorizedArrayType>
+          phi(matrix_free,
+              range,
+              true,
+              dof_no,
+              quad_no,
+              first_selected_component);
+        helper.initialize(phi);
+
+        for (unsigned int face = range.first; face < range.second; ++face)
+          {
+            helper.reinit(face);
+
+            for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
+              {
+                helper.prepare_basis_vector(i);
+                boundary_operation(phi);
+                helper.submit();
+              }
+
+            helper.distribute_local_to_global(diagonal_global_components);
+          }
+      };
+
+    if (face_operation || boundary_operation)
+      matrix_free.template loop<VectorType, int>(cell_operation_wrapped,
+                                                 face_operation_wrapped,
+                                                 boundary_operation_wrapped,
+                                                 diagonal_global,
+                                                 dummy,
+                                                 false);
+    else
+      matrix_free.template cell_loop<VectorType, int>(cell_operation_wrapped,
+                                                      diagonal_global,
+                                                      dummy,
+                                                      false);
   }
 
   template <typename CLASS,
@@ -1084,6 +1722,26 @@ namespace MatrixFreeTools
                                                n_components,
                                                Number,
                                                VectorizedArrayType> &) const,
+    void (CLASS::*face_operation)(FEFaceEvaluation<dim,
+                                                   fe_degree,
+                                                   n_q_points_1d,
+                                                   n_components,
+                                                   Number,
+                                                   VectorizedArrayType> &,
+                                  FEFaceEvaluation<dim,
+                                                   fe_degree,
+                                                   n_q_points_1d,
+                                                   n_components,
+                                                   Number,
+                                                   VectorizedArrayType> &)
+      const,
+    void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
+                                                       fe_degree,
+                                                       n_q_points_1d,
+                                                       n_components,
+                                                       Number,
+                                                       VectorizedArrayType> &)
+      const,
     const CLASS       *owning_class,
     const unsigned int dof_no,
     const unsigned int quad_no,
@@ -1098,7 +1756,11 @@ namespace MatrixFreeTools
                      VectorType>(
       matrix_free,
       diagonal_global,
-      [&](auto &feeval) { (owning_class->*cell_operation)(feeval); },
+      [&](auto &phi) { (owning_class->*cell_operation)(phi); },
+      [&](auto &phi_m, auto &phi_p) {
+        (owning_class->*face_operation)(phi_m, phi_p);
+      },
+      [&](auto &phi) { (owning_class->*boundary_operation)(phi); },
       dof_no,
       quad_no,
       first_selected_component);
@@ -1170,99 +1832,435 @@ namespace MatrixFreeTools
                                           n_q_points_1d,
                                           n_components,
                                           Number,
-                                          VectorizedArrayType> &)> &local_vmult,
-    const unsigned int                                              dof_no,
-    const unsigned int                                              quad_no,
+                                          VectorizedArrayType> &)>
+                      &cell_operation,
+    const unsigned int dof_no,
+    const unsigned int quad_no,
     const unsigned int first_selected_component)
   {
-    std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
-      constraints_for_matrix;
-    const AffineConstraints<typename MatrixType::value_type> &constraints =
-      internal::create_new_affine_constraints_if_needed(matrix,
-                                                        constraints_in,
-                                                        constraints_for_matrix);
-
-    matrix_free.template cell_loop<MatrixType, MatrixType>(
-      [&](const auto &, auto &dst, const auto &, const auto range) {
-        FEEvaluation<dim,
-                     fe_degree,
-                     n_q_points_1d,
-                     n_components,
-                     Number,
-                     VectorizedArrayType>
-          integrator(
-            matrix_free, range, dof_no, quad_no, first_selected_component);
-
-        const unsigned int dofs_per_cell = integrator.dofs_per_cell;
+    compute_matrix<dim,
+                   fe_degree,
+                   n_q_points_1d,
+                   n_components,
+                   Number,
+                   VectorizedArrayType,
+                   MatrixType>(matrix_free,
+                               constraints_in,
+                               matrix,
+                               cell_operation,
+                               {},
+                               {},
+                               dof_no,
+                               quad_no,
+                               first_selected_component);
+  }
 
-        std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
-        std::vector<types::global_dof_index> dof_indices_mf(dofs_per_cell);
+  template <typename CLASS,
+            int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename MatrixType>
+  void
+  compute_matrix(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const AffineConstraints<Number>                    &constraints,
+    MatrixType                                         &matrix,
+    void (CLASS::*cell_operation)(FEEvaluation<dim,
+                                               fe_degree,
+                                               n_q_points_1d,
+                                               n_components,
+                                               Number,
+                                               VectorizedArrayType> &) const,
+    const CLASS       *owning_class,
+    const unsigned int dof_no,
+    const unsigned int quad_no,
+    const unsigned int first_selected_component)
+  {
+    compute_matrix<dim,
+                   fe_degree,
+                   n_q_points_1d,
+                   n_components,
+                   Number,
+                   VectorizedArrayType,
+                   MatrixType>(
+      matrix_free,
+      constraints,
+      matrix,
+      [&](auto &phi) { (owning_class->*cell_operation)(phi); },
+      dof_no,
+      quad_no,
+      first_selected_component);
+  }
 
-        std::array<FullMatrix<typename MatrixType::value_type>,
-                   VectorizedArrayType::size()>
-          matrices;
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename MatrixType>
+  void
+  compute_matrix(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const AffineConstraints<Number>                    &constraints_in,
+    MatrixType                                         &matrix,
+    const std::function<void(FEEvaluation<dim,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          n_components,
+                                          Number,
+                                          VectorizedArrayType> &)>
+      &cell_operation,
+    const std::function<void(FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &,
+                             FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &)>
+      &face_operation,
+    const std::function<void(FEFaceEvaluation<dim,
+                                              fe_degree,
+                                              n_q_points_1d,
+                                              n_components,
+                                              Number,
+                                              VectorizedArrayType> &)>
+                      &boundary_operation,
+    const unsigned int dof_no,
+    const unsigned int quad_no,
+    const unsigned int first_selected_component)
+  {
+    using FEEvalType = FEEvaluation<dim,
+                                    fe_degree,
+                                    n_q_points_1d,
+                                    n_components,
+                                    Number,
+                                    VectorizedArrayType>;
+
+    using FEFaceEvalType = FEFaceEvaluation<dim,
+                                            fe_degree,
+                                            n_q_points_1d,
+                                            n_components,
+                                            Number,
+                                            VectorizedArrayType>;
+
+    internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
+      data_cell;
+
+    data_cell.dof_numbers               = {dof_no};
+    data_cell.quad_numbers              = {dof_no};
+    data_cell.first_selected_components = {first_selected_component};
+    data_cell.batch_type                = {0};
+
+    data_cell.op_create =
+      [&](const std::pair<unsigned int, unsigned int> &range) {
+        std::vector<
+          std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, false>>>
+          phi;
+
+        if (!internal::is_fe_nothing<false>(matrix_free,
+                                            range,
+                                            dof_no,
+                                            quad_no,
+                                            first_selected_component,
+                                            fe_degree,
+                                            n_q_points_1d))
+          phi.emplace_back(std::make_unique<FEEvalType>(
+            matrix_free, range, dof_no, quad_no, first_selected_component));
+
+        return phi;
+      };
 
-        std::fill_n(matrices.begin(),
-                    VectorizedArrayType::size(),
-                    FullMatrix<typename MatrixType::value_type>(dofs_per_cell,
-                                                                dofs_per_cell));
+    data_cell.op_reinit = [](auto &phi, const unsigned batch) {
+      static_cast<FEEvalType &>(*phi[0]).reinit(batch);
+    };
 
-        const auto lexicographic_numbering =
-          matrix_free
-            .get_shape_info(dof_no,
-                            quad_no,
-                            first_selected_component,
-                            integrator.get_active_fe_index(),
-                            integrator.get_active_quadrature_index())
-            .lexicographic_numbering;
+    if (cell_operation)
+      data_cell.op_compute = [&](auto &phi) {
+        cell_operation(static_cast<FEEvalType &>(*phi[0]));
+      };
 
-        for (auto cell = range.first; cell < range.second; ++cell)
+    internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+      data_face;
+
+    data_face.dof_numbers               = {dof_no, dof_no};
+    data_face.quad_numbers              = {dof_no, quad_no};
+    data_face.first_selected_components = {first_selected_component,
+                                           first_selected_component};
+    data_face.batch_type                = {1, 2};
+
+    data_face.op_create =
+      [&](const std::pair<unsigned int, unsigned int> &range) {
+        std::vector<
+          std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
+          phi;
+
+        if (!internal::is_fe_nothing<true>(matrix_free,
+                                           range,
+                                           dof_no,
+                                           quad_no,
+                                           first_selected_component,
+                                           fe_degree,
+                                           n_q_points_1d,
+                                           true) &&
+            !internal::is_fe_nothing<true>(matrix_free,
+                                           range,
+                                           dof_no,
+                                           quad_no,
+                                           first_selected_component,
+                                           fe_degree,
+                                           n_q_points_1d,
+                                           false))
           {
-            integrator.reinit(cell);
+            phi.emplace_back(
+              std::make_unique<FEFaceEvalType>(matrix_free,
+                                               range,
+                                               true,
+                                               dof_no,
+                                               quad_no,
+                                               first_selected_component));
+            phi.emplace_back(
+              std::make_unique<FEFaceEvalType>(matrix_free,
+                                               range,
+                                               false,
+                                               dof_no,
+                                               quad_no,
+                                               first_selected_component));
+          }
 
-            const unsigned int n_filled_lanes =
-              matrix_free.n_active_entries_per_cell_batch(cell);
+        return phi;
+      };
 
-            for (unsigned int v = 0; v < n_filled_lanes; ++v)
-              matrices[v] = 0.0;
+    data_face.op_reinit = [](auto &phi, const unsigned batch) {
+      static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
+      static_cast<FEFaceEvalType &>(*phi[1]).reinit(batch);
+    };
 
-            for (unsigned int j = 0; j < dofs_per_cell; ++j)
-              {
-                for (unsigned int i = 0; i < dofs_per_cell; ++i)
-                  integrator.begin_dof_values()[i] =
-                    static_cast<Number>(i == j);
+    if (face_operation)
+      data_face.op_compute = [&](auto &phi) {
+        face_operation(static_cast<FEFaceEvalType &>(*phi[0]),
+                       static_cast<FEFaceEvalType &>(*phi[1]));
+      };
+
+    internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+      data_boundary;
+
+    data_boundary.dof_numbers               = {dof_no};
+    data_boundary.quad_numbers              = {dof_no};
+    data_boundary.first_selected_components = {first_selected_component};
+    data_boundary.batch_type                = {1};
+
+    data_boundary
+      .op_create = [&](const std::pair<unsigned int, unsigned int> &range) {
+      std::vector<
+        std::unique_ptr<FEEvaluationData<dim, VectorizedArrayType, true>>>
+        phi;
+
+      if (!internal::is_fe_nothing<true>(matrix_free,
+                                         range,
+                                         dof_no,
+                                         quad_no,
+                                         first_selected_component,
+                                         fe_degree,
+                                         n_q_points_1d,
+                                         true))
+        phi.emplace_back(std::make_unique<FEFaceEvalType>(
+          matrix_free, range, true, dof_no, quad_no, first_selected_component));
+
+      return phi;
+    };
+
+    data_boundary.op_reinit = [](auto &phi, const unsigned batch) {
+      static_cast<FEFaceEvalType &>(*phi[0]).reinit(batch);
+    };
+
+    if (boundary_operation)
+      data_boundary.op_compute = [&](auto &phi) {
+        boundary_operation(static_cast<FEFaceEvalType &>(*phi[0]));
+      };
+
+    internal::compute_matrix(
+      matrix_free, constraints_in, data_cell, data_face, data_boundary, matrix);
+  }
+
+  namespace internal
+  {
+    template <int dim,
+              typename Number,
+              typename VectorizedArrayType,
+              typename MatrixType>
+    void
+    compute_matrix(
+      const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+      const AffineConstraints<Number>                    &constraints_in,
+      const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, false>
+        &data_cell,
+      const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+        &data_face,
+      const internal::ComputeMatrixScratchData<dim, VectorizedArrayType, true>
+                 &data_boundary,
+      MatrixType &matrix)
+    {
+      std::unique_ptr<AffineConstraints<typename MatrixType::value_type>>
+        constraints_for_matrix;
+      const AffineConstraints<typename MatrixType::value_type> &constraints =
+        internal::create_new_affine_constraints_if_needed(
+          matrix, constraints_in, constraints_for_matrix);
+
+      const auto batch_operation =
+        [&matrix_free, &constraints, &matrix](
+          auto &data, const std::pair<unsigned int, unsigned int> &range) {
+          if (!data.op_compute)
+            return; // nothing to do
+
+          auto phi = data.op_create(range);
+
+          const unsigned int n_blocks = phi.size();
+
+          if (n_blocks == 0)
+            return;
+
+          Table<1, unsigned int> dofs_per_cell(n_blocks);
+
+          Table<1, std::vector<types::global_dof_index>> dof_indices(n_blocks);
+          Table<2, std::vector<types::global_dof_index>> dof_indices_mf(
+            n_blocks, VectorizedArrayType::size());
+          Table<1, std::vector<unsigned int>> lexicographic_numbering(n_blocks);
+          Table<2,
+                std::array<FullMatrix<typename MatrixType::value_type>,
+                           VectorizedArrayType::size()>>
+            matrices(n_blocks, n_blocks);
+
+          for (unsigned int b = 0; b < n_blocks; ++b)
+            {
+              const auto &shape_info = matrix_free.get_shape_info(
+                data.dof_numbers[b],
+                data.quad_numbers[b],
+                data.first_selected_components[b],
+                phi[b]->get_active_fe_index(),
+                phi[b]->get_active_quadrature_index());
+
+              dofs_per_cell[b] =
+                shape_info.dofs_per_component_on_cell * shape_info.n_components;
+
+              dof_indices[b].resize(dofs_per_cell[b]);
+
+              for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+                dof_indices_mf[b][v].resize(dofs_per_cell[b]);
+
+              lexicographic_numbering[b] = shape_info.lexicographic_numbering;
+            }
+
+          for (unsigned int bj = 0; bj < n_blocks; ++bj)
+            for (unsigned int bi = 0; bi < n_blocks; ++bi)
+              std::fill_n(matrices[bi][bj].begin(),
+                          VectorizedArrayType::size(),
+                          FullMatrix<typename MatrixType::value_type>(
+                            dofs_per_cell[bi], dofs_per_cell[bj]));
+
+          for (auto batch = range.first; batch < range.second; ++batch)
+            {
+              data.op_reinit(phi, batch);
+
+              const unsigned int n_filled_lanes =
+                data.is_face ?
+                  matrix_free.n_active_entries_per_face_batch(batch) :
+                  matrix_free.n_active_entries_per_cell_batch(batch);
+
+              for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                for (unsigned int b = 0; b < n_blocks; ++b)
+                  {
+                    unsigned int const cell_index =
+                      (data.batch_type[b] == 0) ?
+                        (batch * VectorizedArrayType::size() + v) :
+                        ((data.batch_type[b] == 1) ?
+                           matrix_free.get_face_info(batch).cells_interior[v] :
+                           matrix_free.get_face_info(batch).cells_exterior[v]);
+
+                    const auto cell_iterator = matrix_free.get_cell_iterator(
+                      cell_index / VectorizedArrayType::size(),
+                      cell_index % VectorizedArrayType::size(),
+                      data.dof_numbers[b]);
+
+                    if (matrix_free.get_mg_level() !=
+                        numbers::invalid_unsigned_int)
+                      cell_iterator->get_mg_dof_indices(dof_indices[b]);
+                    else
+                      cell_iterator->get_dof_indices(dof_indices[b]);
+
+                    for (unsigned int j = 0; j < dof_indices[b].size(); ++j)
+                      dof_indices_mf[b][v][j] =
+                        dof_indices[b][lexicographic_numbering[b][j]];
+                  }
 
-                local_vmult(integrator);
+              for (unsigned int bj = 0; bj < n_blocks; ++bj)
+                {
+                  for (unsigned int j = 0; j < dofs_per_cell[bj]; ++j)
+                    {
+                      for (unsigned int bi = 0; bi < n_blocks; ++bi)
+                        for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
+                          phi[bi]->begin_dof_values()[i] =
+                            (bj == bi) ? static_cast<Number>(i == j) : 0.0;
+
+                      data.op_compute(phi);
+
+                      for (unsigned int bi = 0; bi < n_blocks; ++bi)
+                        for (unsigned int i = 0; i < dofs_per_cell[bi]; ++i)
+                          for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                            matrices[bi][bj][v](i, j) =
+                              phi[bi]->begin_dof_values()[i][v];
+                    }
 
-                for (unsigned int i = 0; i < dofs_per_cell; ++i)
                   for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                    matrices[v](i, j) = integrator.begin_dof_values()[i][v];
-              }
+                    for (unsigned int bi = 0; bi < n_blocks; ++bi)
+                      constraints.distribute_local_to_global(
+                        matrices[bi][bj][v],
+                        dof_indices_mf[bi][v],
+                        dof_indices_mf[bj][v],
+                        matrix);
+                }
+            }
+        };
 
-            for (unsigned int v = 0; v < n_filled_lanes; ++v)
-              {
-                const auto cell_v =
-                  matrix_free.get_cell_iterator(cell, v, dof_no);
+      const auto cell_operation_wrapped =
+        [&](const auto &, auto &, const auto &, const auto range) {
+          batch_operation(data_cell, range);
+        };
 
-                if (matrix_free.get_mg_level() != numbers::invalid_unsigned_int)
-                  cell_v->get_mg_dof_indices(dof_indices);
-                else
-                  cell_v->get_dof_indices(dof_indices);
+      const auto face_operation_wrapped =
+        [&](const auto &, auto &, const auto &, const auto range) {
+          batch_operation(data_face, range);
+        };
 
-                for (unsigned int j = 0; j < dof_indices.size(); ++j)
-                  dof_indices_mf[j] = dof_indices[lexicographic_numbering[j]];
+      const auto boundary_operation_wrapped =
+        [&](const auto &, auto &, const auto &, const auto range) {
+          batch_operation(data_boundary, range);
+        };
 
-                constraints.distribute_local_to_global(matrices[v],
-                                                       dof_indices_mf,
-                                                       dst);
-              }
-          }
-      },
-      matrix,
-      matrix);
+      if (data_face.op_compute || data_boundary.op_compute)
+        {
+          matrix_free.template loop<MatrixType, MatrixType>(
+            cell_operation_wrapped,
+            face_operation_wrapped,
+            boundary_operation_wrapped,
+            matrix,
+            matrix);
+        }
+      else
+        matrix_free.template cell_loop<MatrixType, MatrixType>(
+          cell_operation_wrapped, matrix, matrix);
 
-    matrix.compress(VectorOperation::add);
-  }
+      matrix.compress(VectorOperation::add);
+    }
+  } // namespace internal
 
   template <typename CLASS,
             int dim,
@@ -1283,6 +2281,26 @@ namespace MatrixFreeTools
                                                n_components,
                                                Number,
                                                VectorizedArrayType> &) const,
+    void (CLASS::*face_operation)(FEFaceEvaluation<dim,
+                                                   fe_degree,
+                                                   n_q_points_1d,
+                                                   n_components,
+                                                   Number,
+                                                   VectorizedArrayType> &,
+                                  FEFaceEvaluation<dim,
+                                                   fe_degree,
+                                                   n_q_points_1d,
+                                                   n_components,
+                                                   Number,
+                                                   VectorizedArrayType> &)
+      const,
+    void (CLASS::*boundary_operation)(FEFaceEvaluation<dim,
+                                                       fe_degree,
+                                                       n_q_points_1d,
+                                                       n_components,
+                                                       Number,
+                                                       VectorizedArrayType> &)
+      const,
     const CLASS       *owning_class,
     const unsigned int dof_no,
     const unsigned int quad_no,
@@ -1298,7 +2316,11 @@ namespace MatrixFreeTools
       matrix_free,
       constraints,
       matrix,
-      [&](auto &feeval) { (owning_class->*cell_operation)(feeval); },
+      [&](auto &phi) { (owning_class->*cell_operation)(phi); },
+      [&](auto &phi_m, auto &phi_p) {
+        (owning_class->*face_operation)(phi_m, phi_p);
+      },
+      [&](auto &phi) { (owning_class->*boundary_operation)(phi); },
       dof_no,
       quad_no,
       first_selected_component);
diff --git a/tests/matrix_free/compute_diagonal_08.cc b/tests/matrix_free/compute_diagonal_08.cc
new file mode 100644 (file)
index 0000000..8f2eecb
--- /dev/null
@@ -0,0 +1,276 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2022 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// Compute diagonal and matrix for SIP with MatrixFreeTools.
+
+#include <deal.II/fe/fe_dgq.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 "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test()
+{
+  const int fe_degree       = 3;
+  const int n_points        = fe_degree + 1;
+  const int n_components    = 1;
+  using Number              = double;
+  using VectorizedArrayType = VectorizedArray<Number>;
+  using VectorType          = Vector<Number>;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_ball(tria);
+  tria.refine_global(0);
+
+  const FE_DGQ<dim> fe_q(fe_degree);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe_q);
+
+  AffineConstraints<Number> constraints;
+  constraints.close();
+
+  typename MatrixFree<dim, Number, VectorizedArrayType>::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;
+
+  MappingQ<dim> mapping(1);
+  QGauss<1>     quad(fe_degree + 1);
+
+  MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+  matrix_free.reinit(mapping, dof_handler, constraints, quad, additional_data);
+
+  const auto cell_operation = [&](auto &phi) {
+    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);
+  };
+
+  const auto face_operation = [&](auto &phi_m, auto &phi_p) {
+    phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+    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_p.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+  };
+
+  const auto boundary_operation = [&](auto &phi_m) {
+    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);
+  };
+
+  const auto vmult = [&](auto &dst, const auto &src) {
+    matrix_free.template loop<VectorType, VectorType>(
+      [&](
+        const auto &matrix_free, auto &dst, const auto &src, const auto range) {
+        FEEvaluation<dim,
+                     fe_degree,
+                     n_points,
+                     n_components,
+                     Number,
+                     VectorizedArrayType>
+          phi(matrix_free);
+
+        for (unsigned int cell = range.first; cell < range.second; ++cell)
+          {
+            phi.reinit(cell);
+            phi.read_dof_values(src);
+            cell_operation(phi);
+            phi.set_dof_values(dst);
+          }
+      },
+      [&](
+        const auto &matrix_free, auto &dst, const auto &src, const auto range) {
+        FEFaceEvaluation<dim,
+                         fe_degree,
+                         n_points,
+                         n_components,
+                         Number,
+                         VectorizedArrayType>
+          phi_m(matrix_free, true);
+        FEFaceEvaluation<dim,
+                         fe_degree,
+                         n_points,
+                         n_components,
+                         Number,
+                         VectorizedArrayType>
+          phi_p(matrix_free, false);
+
+        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_p.read_dof_values(src);
+            face_operation(phi_m, phi_p);
+            phi_m.distribute_local_to_global(dst);
+            phi_p.distribute_local_to_global(dst);
+          }
+      },
+      [&](const auto &matrix_free,
+          auto       &dst,
+          const auto &src,
+          const auto  face_range) {
+        FEFaceEvaluation<dim,
+                         fe_degree,
+                         n_points,
+                         n_components,
+                         Number,
+                         VectorizedArrayType>
+          phi_m(matrix_free, true);
+
+        for (unsigned int face = face_range.first; face < face_range.second;
+             face++)
+          {
+            phi_m.reinit(face);
+            phi_m.read_dof_values(src);
+            boundary_operation(phi_m);
+            phi_m.distribute_local_to_global(dst);
+          }
+      },
+      dst,
+      src,
+      true);
+  };
+
+  // Compute matrix and diagonal (brute-force)
+  FullMatrix<Number> full_matrix(dof_handler.n_dofs(), dof_handler.n_dofs());
+  VectorType         diagonal(dof_handler.n_dofs());
+
+  for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
+    {
+      VectorType src(dof_handler.n_dofs());
+      VectorType dst(dof_handler.n_dofs());
+
+      src[i] = 1.0;
+
+      vmult(dst, src);
+
+      diagonal[i] = dst[i];
+      for (unsigned int j = 0; j < dof_handler.n_dofs(); ++j)
+        full_matrix[j][i] = dst[j];
+    }
+
+
+  // Compute matrix via MatrixFreeTools
+  FullMatrix<Number> full_matrix_mft(dof_handler.n_dofs(),
+                                     dof_handler.n_dofs());
+
+  MatrixFreeTools::compute_matrix<dim,
+                                  fe_degree,
+                                  n_points,
+                                  n_components,
+                                  Number,
+                                  VectorizedArrayType>(matrix_free,
+                                                       constraints,
+                                                       full_matrix_mft,
+                                                       cell_operation,
+                                                       face_operation,
+                                                       boundary_operation);
+
+  for (unsigned int i = 0; i < full_matrix.m(); ++i)
+    for (unsigned int j = 0; j < full_matrix.n(); ++j)
+      if (std::abs(full_matrix[i][j] - full_matrix_mft[i][j]) > 1e-6)
+        {
+          full_matrix.print_formatted(
+            deallog.get_file_stream(), 3, true, 0, "0.000e+00");
+          deallog << std::endl;
+
+          full_matrix_mft.print_formatted(
+            deallog.get_file_stream(), 3, true, 0, "0.000e+00");
+          deallog << std::endl;
+        }
+
+  // Compute diagonal via MatrixFreeTools
+  VectorType diagonal_mft;
+  matrix_free.initialize_dof_vector(diagonal_mft);
+
+  MatrixFreeTools::compute_diagonal<dim,
+                                    fe_degree,
+                                    n_points,
+                                    n_components,
+                                    Number,
+                                    VectorizedArrayType>(matrix_free,
+                                                         diagonal_mft,
+                                                         cell_operation,
+                                                         face_operation,
+                                                         boundary_operation);
+
+  for (unsigned int i = 0; i < diagonal.size(); ++i)
+    if (std::abs(diagonal[i] - diagonal_mft[i]) > 1e-6)
+      {
+        diagonal.print(deallog.get_file_stream());
+        deallog << std::endl;
+
+        diagonal_mft.print(deallog.get_file_stream());
+        deallog << std::endl;
+      }
+
+  deallog << "OK!" << std::endl;
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  test<2>();
+}
diff --git a/tests/matrix_free/compute_diagonal_08.output b/tests/matrix_free/compute_diagonal_08.output
new file mode 100644 (file)
index 0000000..5055485
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL:0::OK!

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.