]> https://gitweb.dealii.org/ - dealii.git/commitdiff
MGTwoLevelTransfer: use ConstraintInfo also for fine vectors 15158/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 29 Apr 2023 12:39:45 +0000 (14:39 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Sun, 30 Apr 2023 10:33:25 +0000 (12:33 +0200)
include/deal.II/matrix_free/constraint_info.h
include/deal.II/multigrid/mg_transfer_global_coarsening.h
include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h

index 6bf150d581121c33943c67c59a5d725b72a4be27..e2d2d9838daa146ca3ca91b6b8523621a47cc020 100644 (file)
@@ -582,7 +582,8 @@ namespace internal
       const unsigned int n_dofs_per_cell,
       const bool         apply_constraints) const
     {
-      if (apply_constraints == false)
+      if ((row_starts_plain_indices.empty() == false) &&
+          (apply_constraints == false))
         {
           for (unsigned int v = 0; v < n_cells; ++v)
             {
index c03bdddaaec7e8f8c13b5dabba894d64c143d6f6..ce9d9e69bcf45835a08d22fd3603e8591833134b 100644 (file)
@@ -303,8 +303,10 @@ protected:
     const std::shared_ptr<const Utilities::MPI::Partitioner>
       &partitioner_coarse,
     const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner_fine,
-    internal::MatrixFreeFunctions::
-      ConstraintInfo<dim, VectorizedArray<Number, width>> &constraint_info);
+    internal::MatrixFreeFunctions::ConstraintInfo<
+      dim,
+      VectorizedArray<Number, width>> &constraint_info_coarse,
+    std::vector<unsigned int> &        dof_indices_fine);
 
   /**
    * Flag if the finite elements on the fine cells are continuous. If yes,
@@ -361,12 +363,6 @@ protected:
    * are a subset of an external Partitioner object.
    */
   mutable AlignedVector<Number> buffer_fine_embedded;
-
-  /**
-   * DoF indices of the fine cells, expressed in indices local to the MPI
-   * rank.
-   */
-  std::vector<unsigned int> level_dof_indices_fine;
 };
 
 
@@ -429,6 +425,8 @@ template <int dim, typename Number>
 class MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>
   : public MGTwoLevelTransferBase<LinearAlgebra::distributed::Vector<Number>>
 {
+  using VectorizedArrayType = VectorizedArray<Number>;
+
 public:
   /**
    * Set up global coarsening between the given DoFHandler objects (
@@ -584,28 +582,28 @@ private:
     /**
      * Prolongation matrix for non-tensor-product elements.
      */
-    AlignedVector<VectorizedArray<Number>> prolongation_matrix;
+    AlignedVector<VectorizedArrayType> prolongation_matrix;
 
     /**
      * 1d prolongation matrix for tensor-product elements.
      */
-    AlignedVector<VectorizedArray<Number>> prolongation_matrix_1d;
+    AlignedVector<VectorizedArrayType> prolongation_matrix_1d;
 
     /**
      * Restriction matrix for non-tensor-product elements.
      */
-    AlignedVector<VectorizedArray<Number>> restriction_matrix;
+    AlignedVector<VectorizedArrayType> restriction_matrix;
 
     /**
      * 1d restriction matrix for tensor-product elements.
      */
-    AlignedVector<VectorizedArray<Number>> restriction_matrix_1d;
+    AlignedVector<VectorizedArrayType> restriction_matrix_1d;
 
     /**
      * ShapeInfo description of the coarse cell. Needed during the
      * fast application of hanging-node constraints.
      */
-    internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>
+    internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
       shape_info_coarse;
   };
 
@@ -615,22 +613,28 @@ private:
   std::vector<MGTransferScheme> schemes;
 
   /**
-   * Helper class for reading from and writing to global vectors and for
+   * Helper class for reading from and writing to global coarse vectors and for
    * applying constraints.
    */
-  internal::MatrixFreeFunctions::ConstraintInfo<dim, VectorizedArray<Number>>
-    constraint_info;
+  internal::MatrixFreeFunctions::ConstraintInfo<dim, VectorizedArrayType>
+    constraint_info_coarse;
+
+  /**
+   * Helper class for reading from and writing to global fine vectors.
+   */
+  internal::MatrixFreeFunctions::ConstraintInfo<dim, VectorizedArrayType>
+    constraint_info_fine;
 
   /**
    * Weights for continuous elements.
    */
-  std::vector<Number> weights;
+  std::vector<Number> weights; // TODO: vectorize
 
   /**
    * Weights for continuous elements, compressed into 3^dim doubles per
    * cell if possible.
    */
-  AlignedVector<VectorizedArray<Number>> weights_compressed;
+  AlignedVector<VectorizedArrayType> weights_compressed;
 
   /**
    * Number of components.
index 2a01f9a4564138298129aebf02297a2929e8d5bf..d89f92a36b1d90756dab5023cc452eb064d02b3f 100644 (file)
@@ -1115,7 +1115,7 @@ namespace internal
     {
       touch_count.reinit(transfer.partitioner_fine);
 
-      for (const auto i : transfer.level_dof_indices_fine)
+      for (const auto i : transfer.constraint_info_fine.dof_indices)
         touch_count.local_element(i) += 1;
       touch_count.compress(VectorOperation::add);
 
@@ -1455,9 +1455,6 @@ namespace internal
 
       // indices
       {
-        transfer.level_dof_indices_fine.resize(n_dof_indices_fine.back(),
-                                               numbers::invalid_unsigned_int);
-
         std::vector<types::global_dof_index> local_dof_indices(
           transfer.schemes[0].n_dofs_per_cell_coarse);
 
@@ -1486,17 +1483,15 @@ namespace internal
           }
 
         // ------------------------------ indices ------------------------------
-        unsigned int *level_dof_indices_fine_0 =
-          transfer.level_dof_indices_fine.data();
-
-        unsigned int *level_dof_indices_fine_1 =
-          level_dof_indices_fine_0 + transfer.schemes[0].n_dofs_per_cell_fine *
-                                       transfer.schemes[0].n_coarse_cells;
+        std::vector<types::global_dof_index> level_dof_indices_fine_0(
+          transfer.schemes[0].n_dofs_per_cell_fine);
+        std::vector<types::global_dof_index> level_dof_indices_fine_1(
+          transfer.schemes[1].n_dofs_per_cell_fine);
 
         unsigned int cell_no_0 = 0;
         unsigned int cell_no_1 = transfer.schemes[0].n_coarse_cells;
 
-        transfer.constraint_info.reinit(
+        transfer.constraint_info_coarse.reinit(
           dof_handler_coarse,
           transfer.schemes[0].n_coarse_cells +
             transfer.schemes[1].n_coarse_cells,
@@ -1504,12 +1499,16 @@ namespace internal
             use_fast_hanging_node_algorithm(dof_handler_coarse,
                                             mg_level_coarse));
 
+        transfer.constraint_info_fine.reinit(
+          transfer.schemes[0].n_coarse_cells +
+          transfer.schemes[1].n_coarse_cells);
+
         process_cells(
           [&](const auto &cell_coarse, const auto &cell_fine) {
             // parent
             {
-              transfer.constraint_info.read_dof_indices(
-                cell_no_0++,
+              transfer.constraint_info_coarse.read_dof_indices(
+                cell_no_0,
                 mg_level_coarse,
                 cell_coarse,
                 constraints_coarse,
@@ -1523,26 +1522,30 @@ namespace internal
                    i < transfer.schemes[0].n_dofs_per_cell_coarse;
                    i++)
                 level_dof_indices_fine_0[i] =
-                  transfer.partitioner_fine->global_to_local(
-                    local_dof_indices[lexicographic_numbering_fine[i]]);
+                  local_dof_indices[lexicographic_numbering_fine[i]];
+
+              transfer.constraint_info_fine.read_dof_indices(
+                cell_no_0, level_dof_indices_fine_0, transfer.partitioner_fine);
             }
 
             // move pointers
             {
-              level_dof_indices_fine_0 +=
-                transfer.schemes[0].n_dofs_per_cell_fine;
+              cell_no_0++;
             }
           },
           [&](const auto &cell_coarse, const auto &cell_fine, const auto c) {
             // parent (only once at the beginning)
             if (c == 0)
               {
-                transfer.constraint_info.read_dof_indices(
-                  cell_no_1++,
+                transfer.constraint_info_coarse.read_dof_indices(
+                  cell_no_1,
                   mg_level_coarse,
                   cell_coarse,
                   constraints_coarse,
                   transfer.partitioner_coarse);
+
+                level_dof_indices_fine_1.assign(level_dof_indices_fine_1.size(),
+                                                numbers::invalid_dof_index);
               }
 
             // child
@@ -1552,12 +1555,12 @@ namespace internal
                    i < transfer.schemes[1].n_dofs_per_cell_coarse;
                    ++i)
                 {
-                  const auto index = transfer.partitioner_fine->global_to_local(
-                    local_dof_indices[lexicographic_numbering_fine[i]]);
+                  const auto index =
+                    local_dof_indices[lexicographic_numbering_fine[i]];
 
                   Assert(level_dof_indices_fine_1
                                [cell_local_children_indices[c][i]] ==
-                             numbers::invalid_unsigned_int ||
+                             numbers::invalid_dof_index ||
                            level_dof_indices_fine_1
                                [cell_local_children_indices[c][i]] == index,
                          ExcInternalError());
@@ -1570,12 +1573,17 @@ namespace internal
             // move pointers (only once at the end)
             if (c + 1 == GeometryInfo<dim>::max_children_per_cell)
               {
-                level_dof_indices_fine_1 +=
-                  transfer.schemes[1].n_dofs_per_cell_fine;
+                transfer.constraint_info_fine.read_dof_indices(
+                  cell_no_1,
+                  level_dof_indices_fine_1,
+                  transfer.partitioner_fine);
+
+                cell_no_1++;
               }
           });
 
-        transfer.constraint_info.finalize();
+        transfer.constraint_info_coarse.finalize();
+        transfer.constraint_info_fine.finalize();
       }
 
 
@@ -1706,9 +1714,11 @@ namespace internal
           Number *weights_0 = transfer.weights.data() + n_dof_indices_fine[0];
           Number *weights_1 = transfer.weights.data() + n_dof_indices_fine[1];
           unsigned int *dof_indices_fine_0 =
-            transfer.level_dof_indices_fine.data() + n_dof_indices_fine[0];
+            transfer.constraint_info_fine.dof_indices.data() +
+            n_dof_indices_fine[0];
           unsigned int *dof_indices_fine_1 =
-            transfer.level_dof_indices_fine.data() + n_dof_indices_fine[1];
+            transfer.constraint_info_fine.dof_indices.data() +
+            n_dof_indices_fine[1];
 
           process_cells(
             [&](const auto &, const auto &) {
@@ -1904,6 +1914,8 @@ namespace internal
           local_dof_indices_coarse_lex(fe_index_pairs.size());
         std::vector<std::vector<types::global_dof_index>>
           local_dof_indices_fine(fe_index_pairs.size());
+        std::vector<std::vector<types::global_dof_index>>
+          local_dof_indices_fine_lex(fe_index_pairs.size());
 
         for (const auto &fe_index_pair : fe_index_pairs)
           {
@@ -1913,6 +1925,8 @@ namespace internal
               transfer.schemes[fe_index_pair.second].n_dofs_per_cell_coarse);
             local_dof_indices_fine[fe_index_pair.second].resize(
               transfer.schemes[fe_index_pair.second].n_dofs_per_cell_fine);
+            local_dof_indices_fine_lex[fe_index_pair.second].resize(
+              transfer.schemes[fe_index_pair.second].n_dofs_per_cell_fine);
 
             n_dof_indices_fine[fe_index_pair.second + 1] =
               transfer.schemes[fe_index_pair.second].n_dofs_per_cell_fine *
@@ -1983,25 +1997,17 @@ namespace internal
             cell_no[i + 1] += cell_no[i];
           }
 
-        transfer.level_dof_indices_fine.resize(n_dof_indices_fine.back());
-
         // ------------------------------ indices  -----------------------------
-        std::vector<unsigned int *> level_dof_indices_fine(
-          fe_index_pairs.size());
 
-        for (unsigned int i = 0; i < fe_index_pairs.size(); ++i)
-          {
-            level_dof_indices_fine[i] =
-              transfer.level_dof_indices_fine.data() + n_dof_indices_fine[i];
-          }
-
-        transfer.constraint_info.reinit(
+        transfer.constraint_info_coarse.reinit(
           dof_handler_coarse,
           cell_no.back(),
           constraints_coarse.n_constraints() > 0 &&
             use_fast_hanging_node_algorithm(dof_handler_coarse,
                                             mg_level_coarse));
 
+        transfer.constraint_info_fine.reinit(cell_no.back());
+
         process_cells([&](const auto &cell_coarse, const auto &cell_fine) {
           const auto fe_pair_no =
             fe_index_pairs[std::pair<unsigned int, unsigned int>(
@@ -2009,8 +2015,8 @@ namespace internal
 
           // parent
           {
-            transfer.constraint_info.read_dof_indices(
-              cell_no[fe_pair_no]++,
+            transfer.constraint_info_coarse.read_dof_indices(
+              cell_no[fe_pair_no],
               mg_level_coarse,
               cell_coarse,
               constraints_coarse,
@@ -2020,23 +2026,27 @@ namespace internal
           // child
           {
             cell_fine->get_dof_indices(local_dof_indices_fine[fe_pair_no]);
+
             for (unsigned int i = 0;
                  i < transfer.schemes[fe_pair_no].n_dofs_per_cell_fine;
-                 i++)
-              level_dof_indices_fine[fe_pair_no][i] =
-                transfer.partitioner_fine->global_to_local(
-                  local_dof_indices_fine
-                    [fe_pair_no][lexicographic_numbering_fine[fe_pair_no][i]]);
+                 ++i)
+              local_dof_indices_fine_lex[fe_pair_no][i] = local_dof_indices_fine
+                [fe_pair_no][lexicographic_numbering_fine[fe_pair_no][i]];
+
+            transfer.constraint_info_fine.read_dof_indices(
+              cell_no[fe_pair_no],
+              local_dof_indices_fine_lex[fe_pair_no],
+              transfer.partitioner_fine);
           }
 
           // move pointers
           {
-            level_dof_indices_fine[fe_pair_no] +=
-              transfer.schemes[fe_pair_no].n_dofs_per_cell_fine;
+            cell_no[fe_pair_no]++;
           }
         });
 
-        transfer.constraint_info.finalize();
+        transfer.constraint_info_coarse.finalize();
+        transfer.constraint_info_fine.finalize();
       }
 
       // ------------------------- prolongation matrix -------------------------
@@ -2205,7 +2215,8 @@ namespace internal
           for (unsigned int i = 0; i < fe_index_pairs.size(); ++i)
             {
               level_dof_indices_fine_[i] =
-                transfer.level_dof_indices_fine.data() + n_dof_indices_fine[i];
+                transfer.constraint_info_fine.dof_indices.data() +
+                n_dof_indices_fine[i];
               weights_[i] = transfer.weights.data() + n_dof_indices_fine[i];
             }
 
@@ -2643,16 +2654,13 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
     LinearAlgebra::distributed::Vector<Number> &      dst,
     const LinearAlgebra::distributed::Vector<Number> &src) const
 {
-  using VectorizedArrayType = VectorizedArray<Number>;
-
   const unsigned int n_lanes = VectorizedArrayType::size();
 
   AlignedVector<VectorizedArrayType> evaluation_data_fine;
   AlignedVector<VectorizedArrayType> evaluation_data_coarse;
 
-  const unsigned int *indices_fine = this->level_dof_indices_fine.data();
-  const Number *      weights      = nullptr;
-  const VectorizedArray<Number> *weights_compressed = nullptr;
+  const Number *             weights            = nullptr;
+  const VectorizedArrayType *weights_compressed = nullptr;
 
   if (this->fine_element_is_continuous)
     {
@@ -2691,18 +2699,17 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
 
           // read from src vector (similar to FEEvaluation::read_dof_values())
           internal::VectorReader<Number, VectorizedArrayType> reader;
-          constraint_info.read_write_operation(reader,
-                                               src,
-                                               evaluation_data_coarse.data(),
-                                               cell_counter,
-                                               n_lanes_filled,
-                                               scheme.n_dofs_per_cell_coarse,
-                                               true);
-          constraint_info.apply_hanging_node_constraints(
+          constraint_info_coarse.read_write_operation(
+            reader,
+            src,
+            evaluation_data_coarse.data(),
+            cell_counter,
+            n_lanes_filled,
+            scheme.n_dofs_per_cell_coarse,
+            true);
+          constraint_info_coarse.apply_hanging_node_constraints(
             cell_counter, n_lanes_filled, false, evaluation_data_coarse);
 
-          cell_counter += n_lanes_filled;
-
           // ---------------------------- coarse -------------------------------
           if (needs_interpolation)
             for (int c = n_components - 1; c >= 0; --c)
@@ -2723,36 +2730,40 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
             evaluation_data_fine = evaluation_data_coarse; // TODO
           // ------------------------------ fine -------------------------------
 
-          // weight and write into dst vector
+          // weight
           if (this->fine_element_is_continuous &&
               this->weights_compressed.size() > 0)
             {
               internal::
-                weight_fe_q_dofs_by_entity<dim, -1, VectorizedArray<Number>>(
+                weight_fe_q_dofs_by_entity<dim, -1, VectorizedArrayType>(
                   weights_compressed,
                   n_components,
                   scheme.degree_fine + 1,
                   evaluation_data_fine.begin());
               weights_compressed += Utilities::pow(3, dim);
             }
-
-          for (unsigned int v = 0; v < n_lanes_filled; ++v)
+          else if (this->fine_element_is_continuous)
             {
-              if (this->fine_element_is_continuous &&
-                  this->weights_compressed.size() == 0)
-                for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i)
-                  dst.local_element(indices_fine[i]) +=
-                    evaluation_data_fine[i][v] * weights[i];
-              else
-                for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i)
-                  dst.local_element(indices_fine[i]) +=
-                    evaluation_data_fine[i][v];
+              for (unsigned int v = 0; v < n_lanes_filled; ++v)
+                {
+                  for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i)
+                    evaluation_data_fine[i][v] *= weights[i];
+                  weights += scheme.n_dofs_per_cell_fine;
+                }
+            }
 
-              indices_fine += scheme.n_dofs_per_cell_fine;
+          // add into dst vector
+          internal::VectorDistributorLocalToGlobal<Number, VectorizedArrayType>
+            writer;
+          constraint_info_fine.read_write_operation(writer,
+                                                    dst,
+                                                    evaluation_data_fine.data(),
+                                                    cell_counter,
+                                                    n_lanes_filled,
+                                                    scheme.n_dofs_per_cell_fine,
+                                                    false);
 
-              if (this->fine_element_is_continuous)
-                weights += scheme.n_dofs_per_cell_fine;
-            }
+          cell_counter += n_lanes_filled;
         }
     }
 }
@@ -2814,16 +2825,13 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
     LinearAlgebra::distributed::Vector<Number> &      dst,
     const LinearAlgebra::distributed::Vector<Number> &src) const
 {
-  using VectorizedArrayType = VectorizedArray<Number>;
-
   const unsigned int n_lanes = VectorizedArrayType::size();
 
   AlignedVector<VectorizedArrayType> evaluation_data_fine;
   AlignedVector<VectorizedArrayType> evaluation_data_coarse;
 
-  const unsigned int *indices_fine = this->level_dof_indices_fine.data();
-  const Number *      weights      = nullptr;
-  const VectorizedArray<Number> *weights_compressed = nullptr;
+  const Number *             weights            = nullptr;
+  const VectorizedArrayType *weights_compressed = nullptr;
 
   if (this->fine_element_is_continuous)
     {
@@ -2860,36 +2868,37 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
               (scheme.n_coarse_cells - cell) :
               n_lanes;
 
-          // read from source vector and weight
-          for (unsigned int v = 0; v < n_lanes_filled; ++v)
-            {
-              if (this->fine_element_is_continuous &&
-                  this->weights_compressed.size() == 0)
-                for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i)
-                  evaluation_data_fine[i][v] =
-                    src.local_element(indices_fine[i]) * weights[i];
-              else
-                for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i)
-                  evaluation_data_fine[i][v] =
-                    src.local_element(indices_fine[i]);
-
-              indices_fine += scheme.n_dofs_per_cell_fine;
-
-              if (this->fine_element_is_continuous)
-                weights += scheme.n_dofs_per_cell_fine;
-            }
-
+          // read from source vector
+          internal::VectorReader<Number, VectorizedArrayType> reader;
+          constraint_info_fine.read_write_operation(reader,
+                                                    src,
+                                                    evaluation_data_fine.data(),
+                                                    cell_counter,
+                                                    n_lanes_filled,
+                                                    scheme.n_dofs_per_cell_fine,
+                                                    false);
+
+          // weight
           if (this->fine_element_is_continuous &&
               this->weights_compressed.size() > 0)
             {
               internal::
-                weight_fe_q_dofs_by_entity<dim, -1, VectorizedArray<Number>>(
+                weight_fe_q_dofs_by_entity<dim, -1, VectorizedArrayType>(
                   weights_compressed,
                   n_components,
                   scheme.degree_fine + 1,
                   evaluation_data_fine.begin());
               weights_compressed += Utilities::pow(3, dim);
             }
+          else if (this->fine_element_is_continuous)
+            {
+              for (unsigned int v = 0; v < n_lanes_filled; ++v)
+                {
+                  for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i)
+                    evaluation_data_fine[i][v] *= weights[i];
+                  weights += scheme.n_dofs_per_cell_fine;
+                }
+            }
 
           // ------------------------------ fine -------------------------------
           if (needs_interpolation)
@@ -2915,15 +2924,16 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
           // FEEvaluation::distribute_global_to_local())
           internal::VectorDistributorLocalToGlobal<Number, VectorizedArrayType>
             writer;
-          constraint_info.apply_hanging_node_constraints(
+          constraint_info_coarse.apply_hanging_node_constraints(
             cell_counter, n_lanes_filled, true, evaluation_data_coarse);
-          constraint_info.read_write_operation(writer,
-                                               dst,
-                                               evaluation_data_coarse.data(),
-                                               cell_counter,
-                                               n_lanes_filled,
-                                               scheme.n_dofs_per_cell_coarse,
-                                               true);
+          constraint_info_coarse.read_write_operation(
+            writer,
+            dst,
+            evaluation_data_coarse.data(),
+            cell_counter,
+            n_lanes_filled,
+            scheme.n_dofs_per_cell_coarse,
+            true);
 
           cell_counter += n_lanes_filled;
         }
@@ -2938,8 +2948,6 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
   interpolate(LinearAlgebra::distributed::Vector<Number> &      dst,
               const LinearAlgebra::distributed::Vector<Number> &src) const
 {
-  using VectorizedArrayType = VectorizedArray<Number>;
-
   const unsigned int n_lanes = VectorizedArrayType::size();
 
   const bool use_src_inplace = this->vec_fine.size() == 0;
@@ -2964,8 +2972,6 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
   AlignedVector<VectorizedArrayType> evaluation_data_fine;
   AlignedVector<VectorizedArrayType> evaluation_data_coarse;
 
-  const unsigned int *indices_fine = this->level_dof_indices_fine.data();
-
   unsigned int cell_counter = 0;
 
   for (const auto &scheme : schemes)
@@ -2976,7 +2982,6 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
       if (scheme.n_dofs_per_cell_fine == 0 ||
           scheme.n_dofs_per_cell_coarse == 0)
         {
-          indices_fine += scheme.n_coarse_cells * scheme.n_dofs_per_cell_fine;
           cell_counter += scheme.n_coarse_cells;
 
           continue;
@@ -3005,15 +3010,15 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
               (scheme.n_coarse_cells - cell) :
               n_lanes;
 
-          // read from source vector and weight
-          for (unsigned int v = 0; v < n_lanes_filled; ++v)
-            {
-              for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i)
-                evaluation_data_fine[i][v] =
-                  vec_fine_ptr->local_element(indices_fine[i]);
-
-              indices_fine += scheme.n_dofs_per_cell_fine;
-            }
+          // read from source vector
+          internal::VectorReader<Number, VectorizedArrayType> reader;
+          constraint_info_fine.read_write_operation(reader,
+                                                    *vec_fine_ptr,
+                                                    evaluation_data_fine.data(),
+                                                    cell_counter,
+                                                    n_lanes_filled,
+                                                    scheme.n_dofs_per_cell_fine,
+                                                    false);
 
           // ------------------------------ fine -------------------------------
           if (needs_interpolation)
@@ -3038,13 +3043,14 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
           // write into dst vector (similar to
           // FEEvaluation::set_dof_values_plain())
           internal::VectorSetter<Number, VectorizedArrayType> writer;
-          constraint_info.read_write_operation(writer,
-                                               *vec_coarse_ptr,
-                                               evaluation_data_coarse.data(),
-                                               cell_counter,
-                                               n_lanes_filled,
-                                               scheme.n_dofs_per_cell_coarse,
-                                               false);
+          constraint_info_coarse.read_write_operation(
+            writer,
+            *vec_coarse_ptr,
+            evaluation_data_coarse.data(),
+            cell_counter,
+            n_lanes_filled,
+            scheme.n_dofs_per_cell_coarse,
+            false);
 
           cell_counter += n_lanes_filled;
         }
@@ -3121,8 +3127,10 @@ MGTwoLevelTransferBase<LinearAlgebra::distributed::Vector<Number>>::
       &external_partitioner_coarse,
     const std::shared_ptr<const Utilities::MPI::Partitioner>
       &external_partitioner_fine,
-    internal::MatrixFreeFunctions::
-      ConstraintInfo<dim, VectorizedArray<Number, width>> &constraint_info)
+    internal::MatrixFreeFunctions::ConstraintInfo<
+      dim,
+      VectorizedArray<Number, width>> &constraint_info_coarse,
+    std::vector<unsigned int> &        dof_indices_fine)
 {
   if (this->partitioner_coarse->is_globally_compatible(
         *external_partitioner_coarse))
@@ -3135,11 +3143,11 @@ MGTwoLevelTransferBase<LinearAlgebra::distributed::Vector<Number>>::
     {
       this->vec_coarse.reinit(0);
 
-      for (auto &i : constraint_info.dof_indices)
+      for (auto &i : constraint_info_coarse.dof_indices)
         i = external_partitioner_coarse->global_to_local(
           this->partitioner_coarse->local_to_global(i));
 
-      for (auto &i : constraint_info.plain_dof_indices)
+      for (auto &i : constraint_info_coarse.plain_dof_indices)
         i = external_partitioner_coarse->global_to_local(
           this->partitioner_coarse->local_to_global(i));
 
@@ -3161,7 +3169,7 @@ MGTwoLevelTransferBase<LinearAlgebra::distributed::Vector<Number>>::
     {
       this->vec_fine.reinit(0);
 
-      for (auto &i : this->level_dof_indices_fine)
+      for (auto &i : dof_indices_fine)
         i = external_partitioner_fine->global_to_local(
           this->partitioner_fine->local_to_global(i));
 
@@ -3185,7 +3193,10 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
       &external_partitioner_fine)
 {
   this->internal_enable_inplace_operations_if_possible(
-    external_partitioner_coarse, external_partitioner_fine, constraint_info);
+    external_partitioner_coarse,
+    external_partitioner_fine,
+    constraint_info_coarse,
+    constraint_info_fine.dof_indices);
 }
 
 
@@ -3362,8 +3373,8 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
   size += this->vec_fine.memory_consumption();
   size += this->vec_coarse.memory_consumption();
   size += MemoryConsumption::memory_consumption(weights);
-  size += MemoryConsumption::memory_consumption(this->level_dof_indices_fine);
-  size += constraint_info.memory_consumption();
+  size += constraint_info_coarse.memory_consumption();
+  size += constraint_info_fine.memory_consumption();
 
   return size;
 }

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.