mg_transfer_global_coarsening.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2020 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_mg_transfer_global_coarsening_h
#define dealii_mg_transfer_global_coarsening_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_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_transfer_matrix_free.h>
 
#include <deal.II/non_matching/mapping_info.h>
 
 
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
namespace internal
{
  class MGTwoLevelTransferImplementation;
}
 
namespace RepartitioningPolicyTools
{
  template <int dim, int spacedim>
  class Base;
}
 
template <int dim,
          typename Number,
          typename MemorySpace = ::MemorySpace::Host>
class MGTransferMF;
#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
 
namespace MGTransferGlobalCoarseningTools
{
  enum class PolynomialCoarseningSequenceType
  {
    bisect,
    decrease_by_one,
    go_to_one
  };
 
  unsigned int
  create_next_polynomial_coarsening_degree(
    const unsigned int                      degree,
    const PolynomialCoarseningSequenceType &p_sequence);
 
  std::vector<unsigned int>
  create_polynomial_coarsening_sequence(
    const unsigned int                      max_degree,
    const PolynomialCoarseningSequenceType &p_sequence);
 
  template <int dim, int spacedim>
  std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
  create_geometric_coarsening_sequence(
    const Triangulation<dim, spacedim> &tria);
 
  template <int dim, int spacedim>
  std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
  create_geometric_coarsening_sequence(
    Triangulation<dim, spacedim>                         &tria,
    const RepartitioningPolicyTools::Base<dim, spacedim> &policy,
    const bool preserve_fine_triangulation,
    const bool repartition_fine_triangulation);
 
  template <int dim, int spacedim>
  std::vector<std::shared_ptr<const Triangulation<dim, spacedim>>>
  create_geometric_coarsening_sequence(
    const Triangulation<dim, spacedim>                   &tria,
    const RepartitioningPolicyTools::Base<dim, spacedim> &policy,
    const bool repartition_fine_triangulation = false);
 
} // namespace MGTransferGlobalCoarseningTools
 
 
 
template <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.");
 
  using Number = typename VectorType::value_type;
 
  MGTwoLevelTransferBase();
 
  void
  prolongate_and_add(VectorType &dst, const VectorType &src) const;
 
  void
  restrict_and_add(VectorType &dst, const VectorType &src) const;
 
  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;
 
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 <int dim, 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;
 
public:
  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner_coarse;
 
  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner_fine;
 
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;
};
 
 
 
template <int dim, typename VectorType>
class MGTwoLevelTransfer : public MGTwoLevelTransferBase<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_no_fine,
         const MatrixFree<dim, Number> &matrix_free_coarse,
         const unsigned int             dof_no_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;
 
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 MGTransferMF<dim, Number, typename VectorType::memory_space>;
};
 
 
 
template <int dim, typename VectorType>
class MGTwoLevelTransferNonNested : public MGTwoLevelTransferBase<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;
 
  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 MGTransferMF<dim, Number, typename VectorType::memory_space>;
};
 
 
 
template <int dim, typename Number, typename MemorySpace>
class MGTransferMF : public ::MGLevelGlobalTransfer<
                       LinearAlgebra::distributed::Vector<Number, MemorySpace>>
{
public:
  using VectorType = LinearAlgebra::distributed::Vector<Number, MemorySpace>;
 
  MGTransferMF();
 
  template <typename MGTwoLevelTransferObject>
  MGTransferMF(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);
 
  MGTransferMF(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;
 
  void
  assert_dof_handler(const DoFHandler<dim> &dof_handler_out) const;
 
  MGLevelObject<MGTwoLevelTransfer<dim, VectorType>> internal_transfer;
 
  MGLevelObject<ObserverPointer<MGTwoLevelTransferBase<VectorType>>> transfer;
 
  std::vector<std::shared_ptr<const Utilities::MPI::Partitioner>>
    external_partitioners;
};
 
 
 
template <int dim, typename Number>
class MGTransferBlockMF
  : public MGTransferBlockMatrixFreeBase<
      dim,
      Number,
      MGTransferMF<dim, Number, ::MemorySpace::Host>>
{
public:
  MGTransferBlockMF(const MGTransferMF<dim, Number, ::MemorySpace::Host>
                      &transfer_operator);
 
  MGTransferBlockMF() = default;
 
  MGTransferBlockMF(const MGConstrainedDoFs &mg_constrained_dofs);
 
  MGTransferBlockMF(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 MGTransferMF<dim, Number, ::MemorySpace::Host> &
  get_matrix_free_transfer(const unsigned int b) const override;
 
private:
  std::vector<MGTransferMF<dim, Number, ::MemorySpace::Host>>
    transfer_operators_internal;
 
  std::vector<ObserverPointer<
    const MGTransferMF<dim, Number, ::MemorySpace::Host>>>
    transfer_operators;
};
 
 
 
template <int dim, typename VectorType>
using MGTransferGlobalCoarsening = MGTransferMF<dim,
                                                typename VectorType::value_type,
                                                ::MemorySpace::Host>;
 
template <int dim, typename VectorType>
using MGTransferBlockGlobalCoarsening =
  MGTransferBlockMF<dim, typename VectorType::value_type>;
 
 
#ifndef DOXYGEN
 
/* ----------------------- Inline functions --------------------------------- */
 
 
 
template <int dim, typename Number, typename MemorySpace>
template <typename MGTwoLevelTransferObject>
MGTransferMF<dim, Number, MemorySpace>::MGTransferMF(
  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
MGTransferMF<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
MGTransferMF<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<VectorType> &>(
      static_cast<const MGTwoLevelTransferBase<VectorType> &>(
        Utilities::get_underlying_value(transfer[l])));
}
 
 
 
template <int dim, typename Number, typename MemorySpace>
template <class InVector>
void
MGTransferMF<dim, Number, MemorySpace>::initialize_dof_vector(
  const unsigned int level,
  VectorType        &vec,
  const InVector    &vec_reference,
  const bool         omit_zeroing_entries) 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 (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
MGTransferMF<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);
    }
 
  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
MGTransferMF<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
MGTransferMF<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
MGTransferMF<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 = transfer.min_level();
  const unsigned int max_level = transfer.max_level();
 
  AssertDimension(min_level, dst.min_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]);
        }
    }
}
 
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif