Introduce MGTransferBlockGlobalCoarsening

include/deal.II/base/mg_level_object.h

@@ -97,6 +97,12 @@ class MGLevelObject : public Subscriptor
   const Object &
   operator[](const unsigned int level) const;
 
+  /**
+   * Return object on level max.
+   */
+  const Object &
+  back() const;
+
   /**
    * Delete all previous contents of this object and reset its size according
    * to the values of @p new_minlevel and @p new_maxlevel.
@@ -230,6 +236,14 @@ MGLevelObject<Object>::operator[](const unsigned int i) const
 }
 
 
+template <class Object>
+const Object &
+MGLevelObject<Object>::back() const
+{
+  return this->operator[](this->max_level());
+}
+
+
 template <class Object>
 template <class... Args>
 void

include/deal.II/multigrid/mg_transfer_global_coarsening.h

@@ -28,6 +28,7 @@
 #include <deal.II/matrix_free/shape_info.h>
 
 #include <deal.II/multigrid/mg_base.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -616,6 +617,38 @@ class MGTransferGlobalCoarsening : public dealii::MGTransferBase<VectorType>
 
 
 
+/**
+ * This class works with LinearAlgebra::distributed::BlockVector and
+ * performs exactly the same transfer operations for each block as
+ * MGTransferGlobalCoarsening.
+ */
+template <int dim, typename VectorType>
+class MGTransferBlockGlobalCoarsening
+  : public MGTransferBlockMatrixFreeBase<
+      dim,
+      typename VectorType::value_type,
+      MGTransferGlobalCoarsening<dim, VectorType>>
+{
+public:
+  /**
+   * Constructor.
+   */
+  MGTransferBlockGlobalCoarsening(
+    const MGTransferGlobalCoarsening<dim, VectorType> &transfer_operator);
+
+protected:
+  const MGTransferGlobalCoarsening<dim, VectorType> &
+  get_matrix_free_transfer(const unsigned int b) const override;
+
+private:
+  /**
+   * Non-block version of transfer operation.
+   */
+  const MGTransferGlobalCoarsening<dim, VectorType> &transfer_operator;
+};
+
+
+
 #ifndef DOXYGEN
 
 /* ----------------------- Inline functions --------------------------------- */

include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h

@@ -45,6 +45,7 @@
 
 #include <deal.II/multigrid/mg_tools.h>
 #include <deal.II/multigrid/mg_transfer_global_coarsening.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.templates.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -3163,6 +3164,32 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
   return size;
 }
 
+
+
+template <int dim, typename VectorType>
+MGTransferBlockGlobalCoarsening<dim, VectorType>::
+  MGTransferBlockGlobalCoarsening(
+    const MGTransferGlobalCoarsening<dim, VectorType> &transfer_operator)
+  : MGTransferBlockMatrixFreeBase<dim,
+                                  typename VectorType::value_type,
+                                  MGTransferGlobalCoarsening<dim, VectorType>>(
+      true)
+  , transfer_operator(transfer_operator)
+{}
+
+
+
+template <int dim, typename VectorType>
+const MGTransferGlobalCoarsening<dim, VectorType> &
+MGTransferBlockGlobalCoarsening<dim, VectorType>::get_matrix_free_transfer(
+  const unsigned int b) const
+{
+  (void)b;
+  AssertDimension(b, 0);
+  return transfer_operator;
+}
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif

include/deal.II/multigrid/mg_transfer_matrix_free.h

@@ -321,6 +321,10 @@ class MGTransferBlockMatrixFreeBase
   : public MGTransferBase<LinearAlgebra::distributed::BlockVector<Number>>
 {
 public:
+  MGTransferBlockMatrixFreeBase(const bool same_for_all)
+    : same_for_all(same_for_all)
+  {}
+
   /**
    * Prolongate a vector from level <tt>to_level-1</tt> to level
    * <tt>to_level</tt> using the embedding matrices of the underlying finite
@@ -428,15 +432,17 @@ class MGTransferBlockMatrixFreeBase
 
 protected:
   /**
-   * Non-block matrix-free versions of transfer operation.
+   * Return the right non-block transfer operator. Has to be implemented by
+   * the derived class.
    */
-  std::vector<TransferType> matrix_free_transfer_vector;
+  virtual const TransferType &
+  get_matrix_free_transfer(const unsigned int b) const = 0;
 
   /**
    * A flag to indicate whether the same DoFHandler is used for all
    * the components or if each block has its own DoFHandler.
    */
-  bool same_for_all;
+  const bool same_for_all;
 };
 
 
@@ -520,6 +526,16 @@ class MGTransferBlockMatrixFree
    */
   std::size_t
   memory_consumption() const;
+
+protected:
+  const MGTransferMatrixFree<dim, Number> &
+  get_matrix_free_transfer(const unsigned int b) const override;
+
+private:
+  /**
+   * Non-block matrix-free versions of transfer operation.
+   */
+  std::vector<MGTransferMatrixFree<dim, Number>> matrix_free_transfer_vector;
 };
 
 
@@ -631,7 +647,6 @@ MGTransferBlockMatrixFreeBase<dim, Number, TransferType>::copy_to_mg(
   MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &dst,
   const LinearAlgebra::distributed::BlockVector<Number2> &        src) const
 {
-  AssertDimension(matrix_free_transfer_vector.size(), 1);
   Assert(same_for_all,
          ExcMessage(
            "This object was initialized with support for usage with one "
@@ -662,78 +677,26 @@ MGTransferBlockMatrixFreeBase<dim, Number, TransferType>::copy_to_mg(
   const unsigned int min_level = dst.min_level();
   const unsigned int max_level = dst.max_level();
 
-  // this function is normally called within the Multigrid class with
-  // dst == defect level block vector. At first run this vector is not
-  // initialized. Do this below:
-  {
-    const parallel::TriangulationBase<dim, spacedim> *tria =
-      (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
-        &(dof_handler[0]->get_triangulation())));
-    for (unsigned int i = 1; i < n_blocks; ++i)
-      AssertThrow(
-        (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
-           &(dof_handler[0]->get_triangulation())) == tria),
-        ExcMessage("The DoFHandler use different Triangulations!"));
-
-    MGLevelObject<bool> do_reinit;
-    do_reinit.resize(min_level, max_level);
-    for (unsigned int level = min_level; level <= max_level; ++level)
-      {
-        do_reinit[level] = false;
-        if (dst[level].n_blocks() != n_blocks)
-          {
-            do_reinit[level] = true;
-            continue; // level
-          }
-        for (unsigned int b = 0; b < n_blocks; ++b)
-          {
-            LinearAlgebra::distributed::Vector<Number> &v = dst[level].block(b);
-            if (v.size() !=
-                  dof_handler[b]->locally_owned_mg_dofs(level).size() ||
-                v.locally_owned_size() !=
-                  dof_handler[b]->locally_owned_mg_dofs(level).n_elements())
-              {
-                do_reinit[level] = true;
-                break; // b
-              }
-          }
-      }
-
-    for (unsigned int level = min_level; level <= max_level; ++level)
-      {
-        if (do_reinit[level])
-          {
-            dst[level].reinit(n_blocks);
-            for (unsigned int b = 0; b < n_blocks; ++b)
-              {
-                LinearAlgebra::distributed::Vector<Number> &v =
-                  dst[level].block(b);
-                v.reinit(dof_handler[b]->locally_owned_mg_dofs(level),
-                         dof_handler[b]->get_communicator());
-              }
-            dst[level].collect_sizes();
-          }
-        else
-          dst[level] = 0;
-      }
-  }
+  for (unsigned int level = min_level; level <= max_level; ++level)
+    if (dst[level].n_blocks() != n_blocks)
+      dst[level].reinit(n_blocks);
 
   // FIXME: this a quite ugly as we need a temporary object:
   MGLevelObject<LinearAlgebra::distributed::Vector<Number>> dst_non_block(
     min_level, max_level);
 
   for (unsigned int b = 0; b < n_blocks; ++b)
     {
-      for (unsigned int l = min_level; l <= max_level; ++l)
-        dst_non_block[l].reinit(dst[l].block(b));
       const unsigned int data_block = same_for_all ? 0 : b;
-      matrix_free_transfer_vector[data_block].copy_to_mg(*dof_handler[b],
-                                                         dst_non_block,
-                                                         src.block(b));
+      get_matrix_free_transfer(data_block)
+        .copy_to_mg(*dof_handler[b], dst_non_block, src.block(b));
 
       for (unsigned int l = min_level; l <= max_level; ++l)
         dst[l].block(b) = dst_non_block[l];
     }
+
+  for (unsigned int level = min_level; level <= max_level; ++level)
+    dst[level].collect_sizes();
 }
 
 template <int dim, typename Number, typename TransferType>
@@ -745,7 +708,11 @@ MGTransferBlockMatrixFreeBase<dim, Number, TransferType>::copy_from_mg(
   const MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &src)
   const
 {
-  AssertDimension(matrix_free_transfer_vector.size(), 1);
+  Assert(same_for_all,
+         ExcMessage(
+           "This object was initialized with support for usage with one "
+           "DoFHandler for each block, but this method assumes that "
+           "the same DoFHandler is used for all the blocks!"));
   const std::vector<const DoFHandler<dim, spacedim> *> mg_dofs(dst.n_blocks(),
                                                                &dof_handler);
 
@@ -785,9 +752,8 @@ MGTransferBlockMatrixFreeBase<dim, Number, TransferType>::copy_from_mg(
           src_non_block[l] = src[l].block(b);
         }
       const unsigned int data_block = same_for_all ? 0 : b;
-      matrix_free_transfer_vector[data_block].copy_from_mg(*dof_handler[b],
-                                                           dst.block(b),
-                                                           src_non_block);
+      get_matrix_free_transfer(data_block)
+        .copy_from_mg(*dof_handler[b], dst.block(b), src_non_block);
     }
 }
 

include/deal.II/multigrid/mg_transfer_matrix_free.templates.h

@@ -0,0 +1,86 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_mg_transfer_matrix_free_templates_h
+#define dealii_mg_transfer_matrix_free_templates_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+template <int dim, typename Number, typename TransferType>
+void
+MGTransferBlockMatrixFreeBase<dim, Number, TransferType>::prolongate(
+  const unsigned int                                     to_level,
+  LinearAlgebra::distributed::BlockVector<Number> &      dst,
+  const LinearAlgebra::distributed::BlockVector<Number> &src) const
+{
+  const unsigned int n_blocks = src.n_blocks();
+  AssertDimension(dst.n_blocks(), n_blocks);
+
+  for (unsigned int b = 0; b < n_blocks; ++b)
+    {
+      const unsigned int data_block = same_for_all ? 0 : b;
+      get_matrix_free_transfer(data_block)
+        .prolongate(to_level, dst.block(b), src.block(b));
+    }
+}
+
+
+
+template <int dim, typename Number, typename TransferType>
+void
+MGTransferBlockMatrixFreeBase<dim, Number, TransferType>::prolongate_and_add(
+  const unsigned int                                     to_level,
+  LinearAlgebra::distributed::BlockVector<Number> &      dst,
+  const LinearAlgebra::distributed::BlockVector<Number> &src) const
+{
+  const unsigned int n_blocks = src.n_blocks();
+  AssertDimension(dst.n_blocks(), n_blocks);
+
+  for (unsigned int b = 0; b < n_blocks; ++b)
+    {
+      const unsigned int data_block = same_for_all ? 0 : b;
+      get_matrix_free_transfer(data_block)
+        .prolongate_and_add(to_level, dst.block(b), src.block(b));
+    }
+}
+
+
+
+template <int dim, typename Number, typename TransferType>
+void
+MGTransferBlockMatrixFreeBase<dim, Number, TransferType>::restrict_and_add(
+  const unsigned int                                     from_level,
+  LinearAlgebra::distributed::BlockVector<Number> &      dst,
+  const LinearAlgebra::distributed::BlockVector<Number> &src) const
+{
+  const unsigned int n_blocks = src.n_blocks();
+  AssertDimension(dst.n_blocks(), n_blocks);
+
+  for (unsigned int b = 0; b < n_blocks; ++b)
+    {
+      const unsigned int data_block = same_for_all ? 0 : b;
+      get_matrix_free_transfer(data_block)
+        .restrict_and_add(from_level, dst.block(b), src.block(b));
+    }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

source/multigrid/mg_transfer_global_coarsening.inst.in

@@ -19,6 +19,14 @@ for (deal_II_dimension : DIMENSIONS; S1 : REAL_SCALARS)
   {
     template class MGTwoLevelTransfer<deal_II_dimension,
                                       LinearAlgebra::distributed::Vector<S1>>;
+    template class MGTransferBlockMatrixFreeBase<
+      deal_II_dimension,
+      S1,
+      MGTransferGlobalCoarsening<deal_II_dimension,
+                                 LinearAlgebra::distributed::Vector<S1>>>;
+    template class MGTransferBlockGlobalCoarsening<
+      deal_II_dimension,
+      LinearAlgebra::distributed::Vector<S1>>;
   }
 
 for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)