mg_transfer_matrix_free.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2016 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_mg_transfer_matrix_free_h
#define dealii_mg_transfer_matrix_free_h
 
#include <deal.II/base/memory_space.h>
#include <deal.II/base/mg_level_object.h>
#include <deal.II/base/mpi_remote_point_evaluation.h>
#include <deal.II/base/vectorization.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
 
#include <deal.II/matrix_free/constraint_info.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/shape_info.h>
 
#include <deal.II/multigrid/mg_base.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/multigrid/mg_transfer.h>
#include <deal.II/multigrid/mg_transfer_internal.h>
 
#include <deal.II/non_matching/mapping_info.h>
 
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
namespace internal
{
  class MGTwoLevelTransferImplementation;
} // namespace internal
 
namespace Portable
{
  template <int dim, typename VectorType>
  class MGTwoLevelTransfer;
 
  namespace internal
  {
    class MGTwoLevelTransferImplementation;
  } // namespace internal
} // namespace Portable
 
template <int dim,
          typename Number,
          typename MemorySpace = ::MemorySpace::Host>
class MGTransferMatrixFree;
#endif
 
namespace mg
{
  struct SignalsNonNested
  {
    boost::signals2::signal<void(const bool before)> prolongation_cell_loop;
 
    boost::signals2::signal<void(const bool before)> restriction_cell_loop;
 
    boost::signals2::signal<void(const bool before)> prolongation;
 
    boost::signals2::signal<void(const bool before)> restriction;
  };
} // namespace mg
 
 
 
template <int dim, typename VectorType>
class MGTwoLevelTransferBase : public EnableObserverPointer
{
public:
  static_assert(
    std::is_same_v<
      VectorType,
      LinearAlgebra::distributed::Vector<typename VectorType::value_type,
                                         MemorySpace::Host>> ||
      std::is_same_v<
        VectorType,
        LinearAlgebra::distributed::Vector<typename VectorType::value_type,
                                           MemorySpace::Default>>,
    "This class is currently only implemented for vectors of "
    "type LinearAlgebra::distributed::Vector.");
 
  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner_coarse;
 
  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner_fine;
 
  virtual void
  prolongate_and_add(VectorType &dst, const VectorType &src) const = 0;
 
  virtual void
  restrict_and_add(VectorType &dst, const VectorType &src) const = 0;
 
  virtual void
  interpolate(VectorType &dst, const VectorType &src) const = 0;
 
  virtual std::pair<bool, bool>
  enable_inplace_operations_if_possible(
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_coarse,
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_fine) = 0;
 
  virtual std::size_t
  memory_consumption() const = 0;
 
 
  virtual std::pair<const DoFHandler<dim> *, unsigned int>
  get_dof_handler_fine() const = 0;
};
 
 
 
namespace internal
{
  template <int dim, typename VectorType>
  class MGTwoLevelTransferCore : public MGTwoLevelTransferBase<dim, VectorType>
  {
  public:
    using Number = typename VectorType::value_type;
 
    MGTwoLevelTransferCore();
 
    void
    prolongate_and_add(VectorType &dst, const VectorType &src) const override;
 
    void
    restrict_and_add(VectorType &dst, const VectorType &src) const override;
 
  protected:
    virtual void
    prolongate_and_add_internal(VectorType       &dst,
                                const VectorType &src) const = 0;
 
    virtual void
    restrict_and_add_internal(VectorType &dst, const VectorType &src) const = 0;
 
    void
    update_ghost_values(const VectorType &vec) const;
 
    void
    compress(VectorType &vec, const VectorOperation::values op) const;
 
    void
    zero_out_ghost_values(const VectorType &vec) const;
 
    template <std::size_t width, typename IndexType>
    std::pair<bool, bool>
    internal_enable_inplace_operations_if_possible(
      const std::shared_ptr<const Utilities::MPI::Partitioner>
        &partitioner_coarse,
      const std::shared_ptr<const Utilities::MPI::Partitioner>
           &partitioner_fine,
      bool &vec_fine_needs_ghost_update,
      internal::MatrixFreeFunctions::
        ConstraintInfo<dim, VectorizedArray<Number, width>, IndexType>
                                &constraint_info_coarse,
      std::vector<unsigned int> &dof_indices_fine);
 
    bool fine_element_is_continuous;
 
  protected:
    mutable VectorType vec_coarse;
 
    mutable VectorType vec_fine;
 
    bool vec_fine_needs_ghost_update;
 
    std::shared_ptr<const Utilities::MPI::Partitioner>
      partitioner_coarse_embedded;
 
    std::shared_ptr<const Utilities::MPI::Partitioner>
      partitioner_fine_embedded;
 
    mutable AlignedVector<Number> buffer_coarse_embedded;
 
    mutable AlignedVector<Number> buffer_fine_embedded;
  };
} // namespace internal
 
 
 
template <int dim, typename VectorType>
class MGTwoLevelTransfer
  : public internal::MGTwoLevelTransferCore<dim, VectorType>
{
public:
  static_assert(
    std::is_same_v<
      VectorType,
      LinearAlgebra::distributed::Vector<typename VectorType::value_type,
                                         MemorySpace::Host>> ||
      std::is_same_v<
        VectorType,
        LinearAlgebra::distributed::Vector<typename VectorType::value_type,
                                           MemorySpace::Default>>,
    "This class is currently only implemented for vectors of "
    "type LinearAlgebra::distributed::Vector.");
 
  using Number = typename VectorType::value_type;
 
  using VectorizedArrayType = VectorizedArray<Number>;
 
  void
  reinit_geometric_transfer(
    const DoFHandler<dim>           &dof_handler_fine,
    const DoFHandler<dim>           &dof_handler_coarse,
    const AffineConstraints<Number> &constraint_fine =
      AffineConstraints<Number>(),
    const AffineConstraints<Number> &constraint_coarse =
      AffineConstraints<Number>(),
    const unsigned int mg_level_fine   = numbers::invalid_unsigned_int,
    const unsigned int mg_level_coarse = numbers::invalid_unsigned_int);
 
  void
  reinit_polynomial_transfer(
    const DoFHandler<dim>           &dof_handler_fine,
    const DoFHandler<dim>           &dof_handler_coarse,
    const AffineConstraints<Number> &constraint_fine =
      AffineConstraints<Number>(),
    const AffineConstraints<Number> &constraint_coarse =
      AffineConstraints<Number>(),
    const unsigned int mg_level_fine   = numbers::invalid_unsigned_int,
    const unsigned int mg_level_coarse = numbers::invalid_unsigned_int);
 
  void
  reinit(const DoFHandler<dim>           &dof_handler_fine,
         const DoFHandler<dim>           &dof_handler_coarse,
         const AffineConstraints<Number> &constraint_fine =
           AffineConstraints<Number>(),
         const AffineConstraints<Number> &constraint_coarse =
           AffineConstraints<Number>(),
         const unsigned int mg_level_fine   = numbers::invalid_unsigned_int,
         const unsigned int mg_level_coarse = numbers::invalid_unsigned_int);
 
  void
  reinit(const MatrixFree<dim, Number> &matrix_free_fine,
         const unsigned int             dof_handler_index_fine,
         const MatrixFree<dim, Number> &matrix_free_coarse,
         const unsigned int             dof_handler_index_coarse);
 
  static bool
  fast_polynomial_transfer_supported(const unsigned int fe_degree_fine,
                                     const unsigned int fe_degree_coarse);
 
  void
  interpolate(VectorType &dst, const VectorType &src) const override;
 
  std::pair<bool, bool>
  enable_inplace_operations_if_possible(
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_coarse,
    const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner_fine)
    override;
 
  std::size_t
  memory_consumption() const override;
 
  std::pair<const DoFHandler<dim> *, unsigned int>
  get_dof_handler_fine() const override;
 
protected:
  void
  prolongate_and_add_internal(VectorType       &dst,
                              const VectorType &src) const override;
 
  void
  restrict_and_add_internal(VectorType       &dst,
                            const VectorType &src) const override;
 
private:
  struct MGTransferScheme
  {
    unsigned int n_coarse_cells;
 
    unsigned int n_dofs_per_cell_coarse;
 
    unsigned int n_dofs_per_cell_fine;
 
    unsigned int degree_coarse;
 
    unsigned int degree_fine;
 
    AlignedVector<double> prolongation_matrix;
 
    AlignedVector<double> restriction_matrix;
 
    internal::MatrixFreeFunctions::ShapeInfo<double> shape_info_coarse;
  };
 
  std::vector<MGTransferScheme> schemes;
 
  internal::MatrixFreeFunctions::
    ConstraintInfo<dim, VectorizedArrayType, types::global_dof_index>
      constraint_info_coarse;
 
  internal::MatrixFreeFunctions::
    ConstraintInfo<dim, VectorizedArrayType, types::global_dof_index>
      constraint_info_fine;
 
  struct MatrixFreeRelatedData
  {
    ObserverPointer<const MatrixFree<dim, Number>> matrix_free_fine;
 
    unsigned int dof_handler_index_fine;
 
    ObserverPointer<const MatrixFree<dim, Number>> matrix_free_coarse;
 
    unsigned int dof_handler_index_coarse;
 
    std::vector<std::array<unsigned int, VectorizedArrayType::size()>>
      cell_list_fine_to_coarse;
  };
 
  std::unique_ptr<MatrixFreeRelatedData> matrix_free_data;
 
  std::vector<unsigned int> weights_start;
 
  AlignedVector<VectorizedArrayType> weights;
 
  std::vector<unsigned char> weights_are_compressed;
 
  unsigned int n_components;
 
  ObserverPointer<const DoFHandler<dim>> dof_handler_fine;
 
  unsigned int mg_level_fine;
 
  friend class internal::MGTwoLevelTransferImplementation;
 
  friend class MGTransferMatrixFree<dim, Number, MemorySpace::Host>;
 
  friend class MGTransferMatrixFree<dim, Number, MemorySpace::Default>;
 
  friend class Portable::MGTwoLevelTransfer<
    dim,
    LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>>;
 
  friend class Portable::internal::MGTwoLevelTransferImplementation;
};
 
 
template <int dim, typename VectorType>
class MGTwoLevelTransferCopyToHost
  : public MGTwoLevelTransferBase<dim, VectorType>
{
public:
  static_assert(
    std::is_same_v<
      VectorType,
      LinearAlgebra::distributed::Vector<typename VectorType::value_type,
                                         MemorySpace::Default>>,
    "type LinearAlgebra::distributed::Vector. This class should not be used for host vectors since it would perform redundant transfers. Use MGTwoLevelTransfer instead.");
 
 
  using Number = typename VectorType::value_type;
 
  /*
   * The vector type used by vectors defined on the host.
   */
  using VectorTypeHost =
    LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>;
 
  MGTwoLevelTransferCopyToHost() = default;
 
  void
  reinit(const DoFHandler<dim>           &dof_handler_fine,
         const DoFHandler<dim>           &dof_handler_coarse,
         const AffineConstraints<Number> &constraint_fine =
           AffineConstraints<Number>(),
         const AffineConstraints<Number> &constraint_coarse =
           AffineConstraints<Number>(),
         const unsigned int mg_level_fine   = numbers::invalid_unsigned_int,
         const unsigned int mg_level_coarse = numbers::invalid_unsigned_int);
 
  void
  reinit(const MatrixFree<dim, Number> &matrix_free_fine,
         const unsigned int             dof_handler_index_fine,
         const MatrixFree<dim, Number> &matrix_free_coarse,
         const unsigned int             dof_handler_index_coarse);
 
  bool
  fast_polynomial_transfer_supported(const unsigned int fe_degree_fine,
                                     const unsigned int fe_degree_coarse);
 
  void
  prolongate_and_add(VectorType &dst, const VectorType &src) const override;
 
  void
  restrict_and_add(VectorType &dst, const VectorType &src) const override;
 
  void
  interpolate(VectorType &dst, const VectorType &src) const override;
 
  std::pair<bool, bool>
  enable_inplace_operations_if_possible(
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_coarse,
    const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner_fine)
    override;
 
  std::size_t
  memory_consumption() const override;
 
  std::pair<const DoFHandler<dim> *, unsigned int>
  get_dof_handler_fine() const override;
 
private:
  MGTwoLevelTransfer<dim, VectorTypeHost> host_transfer;
 
  mutable VectorTypeHost host_vector_coarse;
 
  mutable VectorTypeHost host_vector_fine;
 
  void
  copy_to_host(VectorTypeHost &dst, const VectorType &src) const;
 
  void
  copy_from_host(VectorType &dst, const VectorTypeHost &src) const;
 
  friend class internal::MGTwoLevelTransferImplementation;
 
  friend class MGTransferMatrixFree<dim,
                                    Number,
                                    typename VectorType::memory_space>;
};
 
 
 
template <int dim, typename VectorType>
class MGTwoLevelTransferNonNested
  : public internal::MGTwoLevelTransferCore<dim, VectorType>
{
private:
  static_assert(
    std::is_same_v<
      VectorType,
      LinearAlgebra::distributed::Vector<typename VectorType::value_type>>,
    "This class is currently only implemented for vectors of "
    "type LinearAlgebra::distributed::Vector.");
 
  using Number              = typename VectorType::value_type;
  using VectorizedArrayType = VectorizedArray<Number, 1>;
 
  mg::SignalsNonNested signals_non_nested;
 
public:
  struct AdditionalData
  {
    AdditionalData(const double       tolerance                = 1e-6,
                   const unsigned int rtree_level              = 0,
                   const bool         enforce_all_points_found = true)
      : tolerance(tolerance)
      , rtree_level(rtree_level)
      , enforce_all_points_found(enforce_all_points_found)
    {}
 
    double tolerance;
 
    unsigned int rtree_level;
 
    bool enforce_all_points_found;
  };
 
  MGTwoLevelTransferNonNested(const AdditionalData &data = AdditionalData());
 
  void
  reinit(const DoFHandler<dim>           &dof_handler_fine,
         const DoFHandler<dim>           &dof_handler_coarse,
         const Mapping<dim>              &mapping_fine,
         const Mapping<dim>              &mapping_coarse,
         const AffineConstraints<Number> &constraint_fine =
           AffineConstraints<Number>(),
         const AffineConstraints<Number> &constraint_coarse =
           AffineConstraints<Number>());
 
  void
  interpolate(VectorType &dst, const VectorType &src) const override;
 
  std::pair<bool, bool>
  enable_inplace_operations_if_possible(
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_coarse,
    const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner_fine)
    override;
 
  std::size_t
  memory_consumption() const override;
 
  std::pair<const DoFHandler<dim> *, unsigned int>
  get_dof_handler_fine() const override;
 
  boost::signals2::connection
  connect_prolongation_cell_loop(const std::function<void(const bool)> &slot);
 
  boost::signals2::connection
  connect_restriction_cell_loop(const std::function<void(const bool)> &slot);
 
  boost::signals2::connection
  connect_prolongation(const std::function<void(const bool)> &slot);
 
  boost::signals2::connection
  connect_restriction(const std::function<void(const bool)> &slot);
 
protected:
  AdditionalData additional_data;
  void
  prolongate_and_add_internal(VectorType       &dst,
                              const VectorType &src) const override;
 
  void
  restrict_and_add_internal(VectorType       &dst,
                            const VectorType &src) const override;
 
private:
  template <int n_components>
  void
  prolongate_and_add_internal_comp(VectorType       &dst,
                                   const VectorType &src) const;
 
  template <int n_components>
  void
  restrict_and_add_internal_comp(VectorType &dst, const VectorType &src) const;
 
  ObserverPointer<const DoFHandler<dim>> dof_handler_fine;
 
  unsigned int mg_level_fine;
 
  Utilities::MPI::RemotePointEvaluation<dim> rpe;
 
  mutable std::vector<Number> rpe_input_output;
 
  mutable std::vector<Number> rpe_buffer;
 
  std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>> mapping_info;
 
  internal::MatrixFreeFunctions::
    ConstraintInfo<dim, VectorizedArrayType, unsigned int>
      constraint_info;
 
  std::unique_ptr<FiniteElement<dim>> fe_coarse;
 
  std::vector<unsigned int> level_dof_indices_fine;
 
  std::vector<unsigned int> level_dof_indices_fine_ptrs;
 
  friend class MGTransferMatrixFree<dim,
                                    Number,
                                    typename VectorType::memory_space>;
};
 
 
 
template <int dim, typename Number, typename MemorySpace>
class MGTransferMatrixFree
  : public ::MGLevelGlobalTransfer<
      LinearAlgebra::distributed::Vector<Number, MemorySpace>>
{
public:
  using VectorType = LinearAlgebra::distributed::Vector<Number, MemorySpace>;
 
  MGTransferMatrixFree();
 
  template <typename MGTwoLevelTransferObject>
  MGTransferMatrixFree(
    const MGLevelObject<MGTwoLevelTransferObject> &transfer,
    const std::function<void(const unsigned int, VectorType &)>
      &initialize_dof_vector = {});
 
  template <typename MGTwoLevelTransferObject>
  void
  initialize_two_level_transfers(
    const MGLevelObject<MGTwoLevelTransferObject> &transfer);
 
  void
  build(const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
          &external_partitioners = {});
 
  void
  build(const std::function<void(const unsigned int, VectorType &)>
          &initialize_dof_vector);
 
  MGTransferMatrixFree(const MGConstrainedDoFs &mg_constrained_dofs);
 
  void
  initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);
 
  void
  build(const DoFHandler<dim> &dof_handler,
        const std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
          &external_partitioners = {});
 
  void
  build(const DoFHandler<dim> &dof_handler,
        const std::function<void(const unsigned int, VectorType &)>
          &initialize_dof_vector);
 
  void
  prolongate(const unsigned int to_level,
             VectorType        &dst,
             const VectorType  &src) const override;
 
  void
  prolongate_and_add(const unsigned int to_level,
                     VectorType        &dst,
                     const VectorType  &src) const override;
 
  virtual void
  restrict_and_add(const unsigned int from_level,
                   VectorType        &dst,
                   const VectorType  &src) const override;
 
  template <class InVector>
  void
  copy_to_mg(const DoFHandler<dim>     &dof_handler,
             MGLevelObject<VectorType> &dst,
             const InVector            &src) const;
 
  template <class OutVector>
  void
  copy_from_mg(const DoFHandler<dim>           &dof_handler,
               OutVector                       &dst,
               const MGLevelObject<VectorType> &src) const;
 
  template <class InVector>
  void
  interpolate_to_mg(const DoFHandler<dim>     &dof_handler,
                    MGLevelObject<VectorType> &dst,
                    const InVector            &src) const;
 
  template <class InVector>
  void
  interpolate_to_mg(MGLevelObject<VectorType> &dst, const InVector &src) const;
 
  std::size_t
  memory_consumption() const;
 
  unsigned int
  min_level() const;
 
  unsigned int
  max_level() const;
 
  void
  clear();
 
private:
  void
  initialize_internal_transfer(
    const DoFHandler<dim>                          &dof_handler,
    const ObserverPointer<const MGConstrainedDoFs> &mg_constrained_dofs);
 
  std::pair<const DoFHandler<dim> *, unsigned int>
  get_dof_handler_fine() const;
 
  void
  fill_and_communicate_copy_indices_global_coarsening(
    const DoFHandler<dim> &dof_handler);
 
  template <typename MGTwoLevelTransferObject>
  void
  initialize_transfer_references(
    const MGLevelObject<MGTwoLevelTransferObject> &transfer);
 
  template <class InVector>
  void
  initialize_dof_vector(const unsigned int level,
                        VectorType        &vector,
                        const InVector    &vector_reference,
                        const bool         omit_zeroing_entries,
                        const bool force_same_partitioner = false) const;
 
  void
  assert_dof_handler(const DoFHandler<dim> &dof_handler_out) const;
 
  MGLevelObject<MGTwoLevelTransfer<
    dim,
    LinearAlgebra::distributed::Vector<Number, ::MemorySpace::Host>>>
    internal_transfer;
 
  MGLevelObject<ObserverPointer<MGTwoLevelTransferBase<dim, VectorType>>>
    transfer;
 
  std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
    external_partitioners;
};
 
template <int dim, typename Number, typename TransferType>
class MGTransferBlockMatrixFreeBase
  : public MGTransferBase<LinearAlgebra::distributed::BlockVector<Number>>
{
public:
  explicit MGTransferBlockMatrixFreeBase(const bool same_for_all)
    : same_for_all(same_for_all)
  {}
 
  virtual void
  prolongate(
    const unsigned int                                     to_level,
    LinearAlgebra::distributed::BlockVector<Number>       &dst,
    const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
 
  virtual void
  prolongate_and_add(
    const unsigned int                                     to_level,
    LinearAlgebra::distributed::BlockVector<Number>       &dst,
    const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
 
  virtual void
  restrict_and_add(
    const unsigned int                                     from_level,
    LinearAlgebra::distributed::BlockVector<Number>       &dst,
    const LinearAlgebra::distributed::BlockVector<Number> &src) const override;
 
  template <typename BlockVectorType2>
  void
  copy_to_mg(
    const DoFHandler<dim>                                          &dof_handler,
    MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &dst,
    const BlockVectorType2                                         &src) const;
 
  template <typename BlockVectorType2>
  void
  copy_to_mg(
    const std::vector<const DoFHandler<dim> *>                     &dof_handler,
    MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &dst,
    const BlockVectorType2                                         &src) const;
 
  template <typename BlockVectorType2>
  void
  copy_from_mg(
    const DoFHandler<dim> &dof_handler,
    BlockVectorType2      &dst,
    const MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &src)
    const;
 
  template <typename BlockVectorType2>
  void
  copy_from_mg(
    const std::vector<const DoFHandler<dim> *> &dof_handler,
    BlockVectorType2                           &dst,
    const MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &src)
    const;
 
protected:
  virtual const TransferType &
  get_matrix_free_transfer(const unsigned int b) const = 0;
 
  const bool same_for_all;
};
 
template <int dim, typename Number>
class MGTransferBlockMatrixFree
  : public MGTransferBlockMatrixFreeBase<
      dim,
      Number,
      MGTransferMatrixFree<dim, Number, ::MemorySpace::Host>>
{
public:
  MGTransferBlockMatrixFree(
    const MGTransferMatrixFree<dim, Number, ::MemorySpace::Host>
      &transfer_operator);
 
  MGTransferBlockMatrixFree() = default;
 
  MGTransferBlockMatrixFree(const MGConstrainedDoFs &mg_constrained_dofs);
 
  MGTransferBlockMatrixFree(
    const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);
 
  void
  initialize_constraints(const MGConstrainedDoFs &mg_constrained_dofs);
 
  void
  initialize_constraints(
    const std::vector<MGConstrainedDoFs> &mg_constrained_dofs);
 
  void
  build(const DoFHandler<dim> &dof_handler);
 
  void
  build(const std::vector<const DoFHandler<dim> *> &dof_handler);
 
protected:
  const MGTransferMatrixFree<dim, Number, ::MemorySpace::Host> &
  get_matrix_free_transfer(const unsigned int b) const override;
 
private:
  std::vector<MGTransferMatrixFree<dim, Number, ::MemorySpace::Host>>
    transfer_operators_internal;
 
  std::vector<ObserverPointer<
    const MGTransferMatrixFree<dim, Number, ::MemorySpace::Host>>>
    transfer_operators;
};
 
 
#ifndef DOXYGEN
 
/* ----------------------- Inline functions --------------------------------- */
 
template <int dim, typename VectorType>
void
MGTwoLevelTransferCopyToHost<dim, VectorType>::reinit(
  const DoFHandler<dim>           &dof_handler_fine,
  const DoFHandler<dim>           &dof_handler_coarse,
  const AffineConstraints<Number> &constraint_fine,
  const AffineConstraints<Number> &constraint_coarse,
  const unsigned int               mg_level_fine,
  const unsigned int               mg_level_coarse)
{
  host_transfer.reinit(dof_handler_fine,
                       dof_handler_coarse,
                       constraint_fine,
                       constraint_coarse,
                       mg_level_fine,
                       mg_level_coarse);
}
 
template <int dim, typename VectorType>
void
MGTwoLevelTransferCopyToHost<dim, VectorType>::reinit(
  const MatrixFree<dim, Number> &matrix_free_fine,
  const unsigned int             dof_handler_index_fine,
  const MatrixFree<dim, Number> &matrix_free_coarse,
  const unsigned int             dof_handler_index_coarse)
{
  host_transfer.reinit(matrix_free_fine,
                       dof_handler_index_fine,
                       matrix_free_coarse,
                       dof_handler_index_coarse);
}
 
 
template <int dim, typename VectorType>
bool
MGTwoLevelTransferCopyToHost<dim, VectorType>::
  fast_polynomial_transfer_supported(const unsigned int fe_degree_fine,
                                     const unsigned int fe_degree_coarse)
{
  return host_transfer.fast_polynomial_transfer_supported(fe_degree_fine,
                                                          fe_degree_coarse);
}
 
template <int dim, typename VectorType>
void
MGTwoLevelTransferCopyToHost<dim, VectorType>::prolongate_and_add(
  VectorType       &dst,
  const VectorType &src) const
{
  if (host_vector_coarse.size() == 0)
    host_vector_coarse.reinit(src.get_partitioner());
  if (host_vector_fine.size() == 0)
    host_vector_fine.reinit(dst.get_partitioner());
 
  copy_to_host(host_vector_fine, dst);
  copy_to_host(host_vector_coarse, src);
 
  host_transfer.prolongate_and_add(host_vector_fine, host_vector_coarse);
 
  copy_from_host(dst, host_vector_fine);
}
 
template <int dim, typename VectorType>
void
MGTwoLevelTransferCopyToHost<dim, VectorType>::restrict_and_add(
  VectorType       &dst,
  const VectorType &src) const
{
  if (host_vector_coarse.size() == 0)
    host_vector_coarse.reinit(dst.get_partitioner());
  if (host_vector_fine.size() == 0)
    host_vector_fine.reinit(src.get_partitioner());
 
  copy_to_host(host_vector_coarse, dst);
  copy_to_host(host_vector_fine, src);
 
  host_transfer.restrict_and_add(host_vector_coarse, host_vector_fine);
 
  copy_from_host(dst, host_vector_coarse);
}
 
template <int dim, typename VectorType>
void
MGTwoLevelTransferCopyToHost<dim, VectorType>::interpolate(
  VectorType       &dst,
  const VectorType &src) const
{
  if (host_vector_coarse.size() == 0)
    host_vector_coarse.reinit(dst.get_partitioner());
  if (host_vector_fine.size() == 0)
    host_vector_fine.reinit(src.get_partitioner());
 
  copy_to_host(host_vector_fine, src);
 
  host_transfer.interpolate(host_vector_coarse, host_vector_fine);
 
  copy_from_host(dst, host_vector_coarse);
}
 
template <int dim, typename VectorType>
std::pair<bool, bool>
MGTwoLevelTransferCopyToHost<dim, VectorType>::
  enable_inplace_operations_if_possible(
    const std::shared_ptr<const Utilities::MPI::Partitioner>
      &partitioner_coarse,
    const std::shared_ptr<const Utilities::MPI::Partitioner> &partitioner_fine)
{
  return host_transfer.enable_inplace_operations_if_possible(partitioner_coarse,
                                                             partitioner_fine);
}
 
template <int dim, typename VectorType>
std::size_t
MGTwoLevelTransferCopyToHost<dim, VectorType>::memory_consumption() const
{
  return host_transfer.memory_consumption();
}
 
template <int dim, typename VectorType>
std::pair<const DoFHandler<dim> *, unsigned int>
MGTwoLevelTransferCopyToHost<dim, VectorType>::get_dof_handler_fine() const
{
  return host_transfer.get_dof_handler_fine();
}
 
 
template <int dim, typename VectorType>
void
MGTwoLevelTransferCopyToHost<dim, VectorType>::copy_to_host(
  VectorTypeHost   &dst,
  const VectorType &src) const
{
  LinearAlgebra::ReadWriteVector<Number> rw_vector(
    src.get_partitioner()->locally_owned_range());
  rw_vector.import_elements(src, VectorOperation::insert);
 
  dst.import_elements(rw_vector, VectorOperation::insert);
}
 
template <int dim, typename VectorType>
void
MGTwoLevelTransferCopyToHost<dim, VectorType>::copy_from_host(
  VectorType           &dst,
  const VectorTypeHost &src) const
{
  LinearAlgebra::ReadWriteVector<Number> rw_vector(
    src.get_partitioner()->locally_owned_range());
  rw_vector.import_elements(src, VectorOperation::insert);
 
  dst.import_elements(rw_vector, VectorOperation::insert);
}
 
 
 
template <int dim, typename Number, typename MemorySpace>
template <typename MGTwoLevelTransferObject>
MGTransferMatrixFree<dim, Number, MemorySpace>::MGTransferMatrixFree(
  const MGLevelObject<MGTwoLevelTransferObject> &transfer,
  const std::function<void(const unsigned int, VectorType &)>
    &initialize_dof_vector)
{
  this->transfer.clear();
  this->internal_transfer.clear();
 
  this->initialize_transfer_references(transfer);
  this->build(initialize_dof_vector);
}
 
 
 
template <int dim, typename Number, typename MemorySpace>
template <typename MGTwoLevelTransferObject>
void
MGTransferMatrixFree<dim, Number, MemorySpace>::initialize_two_level_transfers(
  const MGLevelObject<MGTwoLevelTransferObject> &transfer)
{
  this->initialize_transfer_references(transfer);
}
 
 
 
template <int dim, typename Number, typename MemorySpace>
template <typename MGTwoLevelTransferObject>
void
MGTransferMatrixFree<dim, Number, MemorySpace>::initialize_transfer_references(
  const MGLevelObject<MGTwoLevelTransferObject> &transfer)
{
  const unsigned int min_level = transfer.min_level();
  const unsigned int max_level = transfer.max_level();
 
  this->transfer.resize(min_level, max_level);
 
  // Note that transfer[min_level] is empty and never used:
  for (unsigned int l = min_level + 1; l <= max_level; ++l)
    this->transfer[l] = &const_cast<MGTwoLevelTransferBase<dim, VectorType> &>(
      static_cast<const MGTwoLevelTransferBase<dim, VectorType> &>(
        Utilities::get_underlying_value(transfer[l])));
}
 
 
 
template <int dim, typename Number, typename MemorySpace>
template <class InVector>
void
MGTransferMatrixFree<dim, Number, MemorySpace>::initialize_dof_vector(
  const unsigned int level,
  VectorType        &vec,
  const InVector    &vec_reference,
  const bool         omit_zeroing_entries,
  const bool         force_same_partitioner) const
{
  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
 
  if (external_partitioners.empty())
    {
      partitioner = vec_reference.get_partitioner();
    }
  else
    {
      Assert(transfer.min_level() <= level && level <= transfer.max_level(),
             ExcInternalError());
 
      partitioner = external_partitioners[level - transfer.min_level()];
    }
 
  // check if vectors are already correctly initialized
 
  // yes: same partitioners are used
  if (vec.get_partitioner().get() == partitioner.get())
    {
      if (omit_zeroing_entries == false)
        vec = 0;
      return; // nothing to do
    }
 
  // yes: vectors are compatible
  if (!force_same_partitioner && vec.size() == partitioner->size() &&
      vec.locally_owned_size() == partitioner->locally_owned_size())
    {
      if (omit_zeroing_entries == false)
        vec = 0;
      return; // nothing to do
    }
 
  // no
  vec.reinit(partitioner, omit_zeroing_entries);
}
 
 
 
template <int dim, typename Number, typename MemorySpace>
template <class InVector>
void
MGTransferMatrixFree<dim, Number, MemorySpace>::copy_to_mg(
  const DoFHandler<dim>     &dof_handler,
  MGLevelObject<VectorType> &dst,
  const InVector            &src) const
{
  assert_dof_handler(dof_handler);
 
  for (unsigned int level = dst.min_level(); level <= dst.max_level(); ++level)
    {
      const bool zero_out_values =
        (this->perform_plain_copy == false &&
         this->perform_renumbered_plain_copy == false) ||
        level != dst.max_level();
 
      this->initialize_dof_vector(level,
                                  dst[level],
                                  src,
                                  !zero_out_values,
                                  /* force_same_partitioner = */ true);
    }
 
  if (this->perform_plain_copy)
    {
      dst[dst.max_level()].copy_locally_owned_data_from(src);
    }
  else if (this->perform_renumbered_plain_copy)
    {
      auto &dst_level = dst[dst.max_level()];
 
      for (unsigned int i = 0; i < this->copy_indices.back().n_cols(); ++i)
        dst_level.local_element(this->copy_indices.back()(1, i)) =
          src.local_element(i);
    }
  else
    {
      this->ghosted_global_vector = src;
      this->ghosted_global_vector.update_ghost_values();
 
      for (unsigned int l = dst.max_level() + 1; l != dst.min_level();)
        {
          --l;
 
          auto &dst_level = dst[l];
 
          const auto copy_unknowns = [&](const auto &indices) {
            for (unsigned int i = 0; i < indices.n_cols(); ++i)
              dst_level.local_element(indices(1, i)) =
                this->ghosted_global_vector.local_element(indices(0, i));
          };
 
          copy_unknowns(this->copy_indices[l]);
          copy_unknowns(this->copy_indices_level_mine[l]);
 
          dst_level.compress(VectorOperation::insert);
        }
    }
}
 
 
 
template <int dim, typename Number, typename MemorySpace>
template <class OutVector>
void
MGTransferMatrixFree<dim, Number, MemorySpace>::copy_from_mg(
  const DoFHandler<dim>           &dof_handler,
  OutVector                       &dst,
  const MGLevelObject<VectorType> &src) const
{
  assert_dof_handler(dof_handler);
 
  if (this->perform_plain_copy)
    {
      dst.zero_out_ghost_values();
      dst.copy_locally_owned_data_from(src[src.max_level()]);
    }
  else if (this->perform_renumbered_plain_copy)
    {
      const auto &src_level = src[src.max_level()];
      dst.zero_out_ghost_values();
      for (unsigned int i = 0; i < this->copy_indices.back().n_cols(); ++i)
        dst.local_element(i) =
          src_level.local_element(this->copy_indices.back()(1, i));
    }
  else
    {
      dst = 0;
      for (unsigned int l = src.min_level(); l <= src.max_level(); ++l)
        {
          auto &ghosted_vector = this->ghosted_level_vector[l];
 
          if (this->ghosted_level_vector[l].size() > 0)
            ghosted_vector = src[l];
 
          const auto *const ghosted_vector_ptr =
            (this->ghosted_level_vector[l].size() > 0) ? &ghosted_vector :
                                                         &src[l];
 
          ghosted_vector_ptr->update_ghost_values();
 
          const auto copy_unknowns = [&](const auto &indices) {
            for (unsigned int i = 0; i < indices.n_cols(); ++i)
              dst.local_element(indices(0, i)) =
                ghosted_vector_ptr->local_element(indices(1, i));
          };
 
          copy_unknowns(this->copy_indices[l]);
          copy_unknowns(this->copy_indices_global_mine[l]);
        }
      dst.compress(VectorOperation::insert);
    }
}
 
 
 
template <int dim, typename Number, typename MemorySpace>
template <class InVector>
void
MGTransferMatrixFree<dim, Number, MemorySpace>::interpolate_to_mg(
  MGLevelObject<VectorType> &dst,
  const InVector            &src) const
{
  DoFHandler<dim> dof_handler_dummy;
 
  this->interpolate_to_mg(dof_handler_dummy, dst, src);
}
 
 
 
template <int dim, typename Number, typename MemorySpace>
template <class InVector>
void
MGTransferMatrixFree<dim, Number, MemorySpace>::interpolate_to_mg(
  const DoFHandler<dim>     &dof_handler,
  MGLevelObject<VectorType> &dst,
  const InVector            &src) const
{
  assert_dof_handler(dof_handler);
 
  const unsigned int min_level = dst.min_level();
  const unsigned int max_level = transfer.max_level();
 
  AssertDimension(max_level, dst.max_level());
 
  for (unsigned int level = min_level; level <= max_level; ++level)
    {
      const bool zero_out_values = false;
      this->initialize_dof_vector(level, dst[level], src, !zero_out_values);
    }
 
  if (this->perform_plain_copy)
    {
      dst[max_level].copy_locally_owned_data_from(src);
 
      for (unsigned int l = max_level; l > min_level; --l)
        this->transfer[l]->interpolate(dst[l - 1], dst[l]);
    }
  else if (this->perform_renumbered_plain_copy)
    {
      auto &dst_level = dst[max_level];
 
      for (unsigned int i = 0; i < this->solution_copy_indices.back().n_cols();
           ++i)
        dst_level.local_element(this->solution_copy_indices.back()(1, i)) =
          src.local_element(i);
 
      for (unsigned int l = max_level; l > min_level; --l)
        this->transfer[l]->interpolate(dst[l - 1], dst[l]);
    }
  else
    {
      this->solution_ghosted_global_vector = src;
      this->solution_ghosted_global_vector.update_ghost_values();
 
      for (unsigned int l = max_level + 1; l != min_level;)
        {
          --l;
 
          auto &dst_level = dst[l];
 
          const auto copy_unknowns = [&](const auto &indices) {
            for (unsigned int i = 0; i < indices.n_cols(); ++i)
              dst_level.local_element(indices(1, i)) =
                this->solution_ghosted_global_vector.local_element(
                  indices(0, i));
          };
 
          copy_unknowns(this->solution_copy_indices[l]);
          copy_unknowns(this->solution_copy_indices_level_mine[l]);
 
          dst_level.compress(VectorOperation::insert);
 
          if (l != min_level)
            this->transfer[l]->interpolate(dst[l - 1], dst[l]);
        }
    }
}
 
 
 
template <int dim, typename Number, typename TransferType>
template <typename BlockVectorType2>
void
MGTransferBlockMatrixFreeBase<dim, Number, TransferType>::copy_to_mg(
  const DoFHandler<dim>                                          &dof_handler,
  MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &dst,
  const BlockVectorType2                                         &src) const
{
  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> *> mg_dofs(src.n_blocks(),
                                                     &dof_handler);
 
  copy_to_mg(mg_dofs, dst, src);
}
 
 
 
template <int dim, typename Number, typename TransferType>
template <typename BlockVectorType2>
void
MGTransferBlockMatrixFreeBase<dim, Number, TransferType>::copy_to_mg(
  const std::vector<const DoFHandler<dim> *>                     &dof_handler,
  MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &dst,
  const BlockVectorType2                                         &src) const
{
  const unsigned int n_blocks = src.n_blocks();
  AssertDimension(dof_handler.size(), n_blocks);
 
  if (n_blocks == 0)
    return;
 
  const unsigned int min_level = dst.min_level();
  const unsigned int max_level = dst.max_level();
 
  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)
    {
      const unsigned int data_block = same_for_all ? 0 : 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>
template <typename BlockVectorType2>
void
MGTransferBlockMatrixFreeBase<dim, Number, TransferType>::copy_from_mg(
  const DoFHandler<dim> &dof_handler,
  BlockVectorType2      &dst,
  const MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &src)
  const
{
  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> *> mg_dofs(dst.n_blocks(),
                                                     &dof_handler);
 
  copy_from_mg(mg_dofs, dst, src);
}
 
template <int dim, typename Number, typename TransferType>
template <typename BlockVectorType2>
void
MGTransferBlockMatrixFreeBase<dim, Number, TransferType>::copy_from_mg(
  const std::vector<const DoFHandler<dim> *> &dof_handler,
  BlockVectorType2                           &dst,
  const MGLevelObject<LinearAlgebra::distributed::BlockVector<Number>> &src)
  const
{
  const unsigned int n_blocks = dst.n_blocks();
  AssertDimension(dof_handler.size(), n_blocks);
 
  if (n_blocks == 0)
    return;
 
  const unsigned int min_level = src.min_level();
  const unsigned int max_level = src.max_level();
 
  for (unsigned int l = min_level; l <= max_level; ++l)
    AssertDimension(src[l].n_blocks(), dst.n_blocks());
 
  // FIXME: this a quite ugly as we need a temporary object:
  MGLevelObject<LinearAlgebra::distributed::Vector<Number>> src_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)
        {
          src_non_block[l].reinit(src[l].block(b));
          src_non_block[l] = src[l].block(b);
        }
      const unsigned int data_block = same_for_all ? 0 : b;
      get_matrix_free_transfer(data_block)
        .copy_from_mg(*dof_handler[b], dst.block(b), src_non_block);
    }
}
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif