Create compute_weights_fe_q_dofs_by_entity()

diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h
index 32cd53edb2605473ec97457e550462d314ed5d5b..fddbbf3ee78b4902f410564d3e8171f7929b7bb8 100644 (file)
--- a/include/deal.II/matrix_free/tensor_product_kernels.h
+++ b/include/deal.II/matrix_free/tensor_product_kernels.h
@@ -3361,6 +3361,67 @@ namespace internal
   }
 
 
+  template <int dim, int n_points_1d_template, typename Number>
+  inline bool
+  compute_weights_fe_q_dofs_by_entity(const Number *     data,
+                                      const unsigned int n_components,
+                                      const int n_points_1d_non_template,
+                                      Number *  weights)
+  {
+    const int n_points_1d = n_points_1d_template != -1 ?
+                              n_points_1d_template :
+                              n_points_1d_non_template;
+
+    Assert(n_points_1d > 0, ExcNotImplemented());
+    Assert(n_points_1d < 100, ExcNotImplemented());
+
+    unsigned int compressed_index[100];
+    compressed_index[0] = 0;
+    for (int i = 1; i < n_points_1d - 1; ++i)
+      compressed_index[i] = 1;
+    compressed_index[n_points_1d - 1] = 2;
+
+    // Insert the number data into a storage position for weight,
+    // ensuring that the array has either not been touched before
+    // or the previous content is the same. In case the previous
+    // content has a different value, we exit this function and
+    // signal to outer functions that the compression was not possible.
+    const auto check_and_set = [](Number &weight, const Number &data) {
+      if (weight == Number(-1.0) || weight == data)
+        {
+          weight = data;
+          return true; // success for the entry
+        }
+
+      return false; // failure for the entry
+    };
+
+    for (unsigned int c = 0; c < Utilities::pow<unsigned int>(3, dim); ++c)
+      weights[c] = Number(-1.0);
+
+    for (unsigned int c = 0; c < n_components; ++c)
+      for (int k = 0; k < (dim > 2 ? n_points_1d : 1); ++k)
+        for (int j = 0; j < (dim > 1 ? n_points_1d : 1);
+             ++j, data += n_points_1d)
+          {
+            const unsigned int shift =
+              9 * compressed_index[k] + 3 * compressed_index[j];
+
+            if (!check_and_set(weights[shift], data[0]))
+              return false; // failure
+
+            for (int i = 1; i < n_points_1d - 1; ++i)
+              if (!check_and_set(weights[shift + 1], data[i]))
+                return false; // failure
+
+            if (!check_and_set(weights[shift + 2], data[n_points_1d - 1]))
+              return false; // failure
+          }
+
+    return true; // success
+  }
+
+
 } // end of namespace internal
 
 

diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
index ce52dda65c15d65b5bf973e59e593e1bdfdf80e0..6fe22ecdcc6c619606284751db0aae9da9ab04f8 100644 (file)
--- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
+++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
@@ -1162,71 +1162,25 @@ namespace internal
       MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>
         &transfer)
     {
-      const Number *weights = transfer.weights.data();
-
-      // Helper function that extracts the weight. Before doing so,
-      //  it checks the values already set.
-      const auto set = [](auto &mask, const auto &weight) {
-        if (mask == -1.0 || mask == weight)
-          {
-            mask = weight;
-            return true;
-          }
-
-        return false;
-      };
-
-      std::vector<std::array<Number, Utilities::pow(3, dim)>> masks;
-
-      // loop over all schemes
+      unsigned int n_cells = 0;
       for (const auto &scheme : transfer.schemes)
-        {
-          const int    loop_length = scheme.degree_fine + 1;
-          unsigned int degree_to_3[100];
-          degree_to_3[0] = 0;
-          for (int i = 1; i < loop_length - 1; ++i)
-            degree_to_3[i] = 1;
-          degree_to_3[loop_length - 1] = 2;
-
-          // loop over all cells
-          for (unsigned int cell = 0; cell < scheme.n_coarse_cells; ++cell)
-            {
-              std::vector<std::array<Number, Utilities::pow(3, dim)>> mask(
-                transfer.n_components);
-
-              // loop over all components
-              for (unsigned int c = 0; c < transfer.n_components; ++c)
-                mask[c].fill(-1.0);
-
-              for (unsigned int c = 0; c < transfer.n_components; ++c)
-                for (int k = 0; k < (dim > 2 ? loop_length : 1); ++k)
-                  for (int j = 0; j < (dim > 1 ? loop_length : 1); ++j)
-                    {
-                      const unsigned int shift =
-                        9 * degree_to_3[k] + 3 * degree_to_3[j];
-
-                      if (!set(mask[c][shift], weights[0]))
-                        return; // early return, since the entry already
-                                // had a different weight
-
-                      for (int i = 1; i < loop_length - 1; ++i)
-                        if (!set(mask[c][shift + 1], weights[i]))
-                          return;
+        n_cells += scheme.n_coarse_cells;
 
-                      if (!set(mask[c][shift + 2], weights[loop_length - 1]))
-                        return;
+      std::vector<std::array<Number, Utilities::pow(3, dim)>> masks(n_cells);
 
-                      weights += loop_length;
-                    }
+      const Number *weight_ptr = transfer.weights.data();
+      std::array<Number, Utilities::pow(3, dim)> *mask_ptr = masks.data();
 
-              if (std::find_if(mask.begin(), mask.end(), [&](const auto &a) {
-                    return a != mask[0];
-                  }) != mask.end())
-                return;
-
-              masks.push_back(mask[0]);
-            }
-        }
+      // (try to) compress weights for each cell
+      for (const auto &scheme : transfer.schemes)
+        for (unsigned int cell = 0; cell < scheme.n_coarse_cells;
+             ++cell, weight_ptr += scheme.n_dofs_per_cell_fine, ++mask_ptr)
+          if (!compute_weights_fe_q_dofs_by_entity<dim, -1, Number>(
+                weight_ptr,
+                transfer.n_components,
+                scheme.degree_fine + 1,
+                mask_ptr->data()))
+            return;
 
       // vectorize weights
       AlignedVector<VectorizedArray<Number>> masks_vectorized;