solution_transfer.cc

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 1999 - 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.
//
// ------------------------------------------------------------------------
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/memory_consumption.h>
 
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/distributed/tria.h>
 
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
 
#include <deal.II/fe/fe.h>
 
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
 
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/petsc_block_vector.h>
#include <deal.II/lac/petsc_vector.h>
#include <deal.II/lac/trilinos_parallel_block_vector.h>
#include <deal.II/lac/trilinos_vector.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_element_access.h>
 
#include <deal.II/numerics/solution_transfer.h>
 
#include <functional>
#include <numeric>
 
DEAL_II_NAMESPACE_OPEN
 
 
namespace
{
  template <typename value_type>
  std::vector<char>
  pack_dof_values(std::vector<Vector<value_type>> &dof_values,
                  const unsigned int               dofs_per_cell)
  {
    for (const auto &values : dof_values)
      {
        AssertDimension(values.size(), dofs_per_cell);
        (void)values;
      }
 
    const std::size_t bytes_per_entry = sizeof(value_type) * dofs_per_cell;
 
    std::vector<char> buffer(dof_values.size() * bytes_per_entry);
    for (unsigned int i = 0; i < dof_values.size(); ++i)
      std::memcpy(&buffer[i * bytes_per_entry],
                  &dof_values[i](0),
                  bytes_per_entry);
 
    return buffer;
  }
 
 
 
  template <typename value_type>
  std::vector<Vector<value_type>>
  unpack_dof_values(
    const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
    const unsigned int dofs_per_cell)
  {
    const std::size_t  bytes_per_entry = sizeof(value_type) * dofs_per_cell;
    const unsigned int n_elements      = data_range.size() / bytes_per_entry;
 
    Assert((data_range.size() % bytes_per_entry == 0), ExcInternalError());
 
    std::vector<Vector<value_type>> unpacked_data;
    unpacked_data.reserve(n_elements);
    for (unsigned int i = 0; i < n_elements; ++i)
      {
        Vector<value_type> dof_values(dofs_per_cell);
        std::memcpy(&dof_values(0),
                    &(*std::next(data_range.begin(), i * bytes_per_entry)),
                    bytes_per_entry);
        unpacked_data.emplace_back(std::move(dof_values));
      }
 
    return unpacked_data;
  }
} // namespace
 
 
 
template <int dim, typename VectorType, int spacedim>
SolutionTransfer<dim, VectorType, spacedim>::SolutionTransfer(
  const DoFHandler<dim, spacedim> &dof,
  const bool                       average_values)
  : dof_handler(&dof, typeid(*this).name())
  , average_values(average_values)
  , handle(numbers::invalid_unsigned_int)
{}
 
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::
  prepare_for_coarsening_and_refinement(
    const std::vector<const VectorType *> &all_in)
{
  const ::internal::parallel::shared::
    TemporarilyRestoreSubdomainIds<dim, spacedim>
      subdomain_modifier(dof_handler->get_triangulation());
 
  for (unsigned int i = 0; i < all_in.size(); ++i)
    Assert(all_in[i]->size() == dof_handler->n_dofs(),
           ExcDimensionMismatch(all_in[i]->size(), dof_handler->n_dofs()));
 
  input_vectors = all_in;
  register_data_attach();
}
 
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::
  prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
{
  std::vector<const VectorType *> temp(all_in.size());
 
  for (std::size_t i = 0; i < temp.size(); ++i)
    temp[i] = &(all_in[i]);
 
  this->prepare_for_coarsening_and_refinement(temp);
}
 
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::register_data_attach()
{
  // TODO: casting away constness is bad
  auto tria = const_cast<::Triangulation<dim, spacedim> *>(
    &dof_handler->get_triangulation());
  Assert(tria != nullptr, ExcInternalError());
 
  Assert(handle == numbers::invalid_unsigned_int,
         ExcMessage("You can only add one solution per "
                    "SolutionTransfer object."));
 
  handle = tria->register_data_attach(
    [this](const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
           const CellStatus                                            status) {
      return this->pack_callback(cell_, status);
    },
    /*returns_variable_size_data=*/dof_handler->has_hp_capabilities());
}
 
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::
  prepare_for_coarsening_and_refinement(const VectorType &in)
{
  std::vector<const VectorType *> all_in(1, &in);
  prepare_for_coarsening_and_refinement(all_in);
}
 
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::prepare_for_serialization(
  const VectorType &in)
{
  std::vector<const VectorType *> all_in(1, &in);
  prepare_for_serialization(all_in);
}
 
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::prepare_for_serialization(
  const std::vector<const VectorType *> &all_in)
{
  prepare_for_coarsening_and_refinement(all_in);
}
 
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::deserialize(VectorType &in)
{
  std::vector<VectorType *> all_in(1, &in);
  deserialize(all_in);
}
 
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::deserialize(
  std::vector<VectorType *> &all_in)
{
  register_data_attach();
 
  // this makes interpolate() happy
  input_vectors.resize(all_in.size());
 
  interpolate(all_in);
}
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::interpolate(
  std::vector<VectorType *> &all_out)
{
  const ::internal::parallel::shared::
    TemporarilyRestoreSubdomainIds<dim, spacedim>
      subdomain_modifier(dof_handler->get_triangulation());
 
  Assert(input_vectors.size() == all_out.size(),
         ExcDimensionMismatch(input_vectors.size(), all_out.size()));
  for (unsigned int i = 0; i < all_out.size(); ++i)
    Assert(all_out[i]->size() == dof_handler->n_dofs(),
           ExcDimensionMismatch(all_out[i]->size(), dof_handler->n_dofs()));
 
  // TODO: casting away constness is bad
  auto tria = const_cast<::Triangulation<dim, spacedim> *>(
    &dof_handler->get_triangulation());
  Assert(tria != nullptr, ExcInternalError());
  Assert(
    handle != numbers::invalid_unsigned_int,
    ExcMessage(
      "You can only call interpolate() once per SolutionTransfer object."));
 
  using Number = typename VectorType::value_type;
 
  if (average_values)
    for (auto *const vec : all_out)
      *vec = Number();
 
  VectorType valence;
 
  // initialize valence vector only if we need to average
  if (average_values)
    valence.reinit(*all_out[0]);
 
  tria->notify_ready_to_unpack(
    handle,
    [this, &all_out, &valence](
      const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
      const CellStatus                                            status,
      const boost::iterator_range<std::vector<char>::const_iterator>
        &data_range) {
      this->unpack_callback(cell_, status, data_range, all_out, valence);
    });
 
  if (average_values)
    {
      // finalize valence: compress and invert
      valence.compress(VectorOperation::add);
      for (const auto i : valence.locally_owned_elements())
        valence[i] = (static_cast<Number>(valence[i]) == Number() ?
                        Number() :
                        (Number(1.0) / static_cast<Number>(valence[i])));
      valence.compress(VectorOperation::insert);
 
      for (auto *const vec : all_out)
        {
          // compress and weight with valence
          vec->compress(VectorOperation::add);
          vec->scale(valence);
        }
    }
  else
    {
      for (auto *const vec : all_out)
        vec->compress(VectorOperation::insert);
    }
 
  input_vectors.clear();
  handle = numbers::invalid_unsigned_int;
}
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::interpolate(
  std::vector<VectorType> &all_out)
{
  std::vector<VectorType *> temp(all_out.size());
 
  for (std::size_t i = 0; i < temp.size(); ++i)
    temp[i] = &(all_out[i]);
 
  this->interpolate(temp);
}
 
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::interpolate(VectorType &out)
{
  std::vector<VectorType *> all_out(1, &out);
  interpolate(all_out);
}
 
 
 
template <int dim, typename VectorType, int spacedim>
std::vector<char>
SolutionTransfer<dim, VectorType, spacedim>::pack_callback(
  const typename Triangulation<dim, spacedim>::cell_iterator &cell_,
  const CellStatus                                            status)
{
  typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_, dof_handler);
 
  // create buffer for each individual object
  std::vector<::Vector<typename VectorType::value_type>> dof_values(
    input_vectors.size());
 
  unsigned int fe_index = 0;
  if (dof_handler->has_hp_capabilities())
    {
      switch (status)
        {
          case CellStatus::cell_will_persist:
          case CellStatus::cell_will_be_refined:
            {
              fe_index = cell->future_fe_index();
              break;
            }
 
          case CellStatus::children_will_be_coarsened:
            {
              // In case of coarsening, we need to find a suitable FE index
              // for the parent cell. We choose the 'least dominant fe'
              // on all children from the associated FECollection.
#ifdef DEBUG
              for (const auto &child : cell->child_iterators())
                Assert(child->is_active() && child->coarsen_flag_set(),
                       typename ::Triangulation<
                         dim>::ExcInconsistentCoarseningFlags());
#endif
 
              fe_index = ::internal::hp::DoFHandlerImplementation::
                dominated_future_fe_on_children<dim, spacedim>(cell);
              break;
            }
 
          default:
            Assert(false, ExcInternalError());
            break;
        }
    }
 
  const unsigned int dofs_per_cell =
    dof_handler->get_fe(fe_index).n_dofs_per_cell();
 
  if (dofs_per_cell == 0)
    return std::vector<char>(); // nothing to do for FE_Nothing
 
  auto it_input  = input_vectors.cbegin();
  auto it_output = dof_values.begin();
  for (; it_input != input_vectors.cend(); ++it_input, ++it_output)
    {
      it_output->reinit(dofs_per_cell);
      cell->get_interpolated_dof_values(*(*it_input), *it_output, fe_index);
    }
 
  return pack_dof_values<typename VectorType::value_type>(dof_values,
                                                          dofs_per_cell);
}
 
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::unpack_callback(
  const typename Triangulation<dim, spacedim>::cell_iterator     &cell_,
  const CellStatus                                                status,
  const boost::iterator_range<std::vector<char>::const_iterator> &data_range,
  std::vector<VectorType *>                                      &all_out,
  VectorType                                                     &valence)
{
  typename DoFHandler<dim, spacedim>::cell_iterator cell(*cell_, dof_handler);
 
  unsigned int fe_index = 0;
  if (dof_handler->has_hp_capabilities())
    {
      switch (status)
        {
          case CellStatus::cell_will_persist:
          case CellStatus::children_will_be_coarsened:
            {
              fe_index = cell->active_fe_index();
              break;
            }
 
          case CellStatus::cell_will_be_refined:
            {
              // After refinement, this particular cell is no longer active,
              // and its children have inherited its FE index. However, to
              // unpack the data on the old cell, we need to recover its FE
              // index from one of the children. Just to be sure, we also
              // check if all children have the same FE index.
              fe_index = cell->child(0)->active_fe_index();
              for (unsigned int child_index = 1;
                   child_index < cell->n_children();
                   ++child_index)
                Assert(cell->child(child_index)->active_fe_index() == fe_index,
                       ExcInternalError());
              break;
            }
 
          default:
            Assert(false, ExcInternalError());
            break;
        }
    }
 
  const unsigned int dofs_per_cell =
    dof_handler->get_fe(fe_index).n_dofs_per_cell();
 
  if (dofs_per_cell == 0)
    return; // nothing to do for FE_Nothing
 
  const std::vector<::Vector<typename VectorType::value_type>>
    dof_values =
      unpack_dof_values<typename VectorType::value_type>(data_range,
                                                         dofs_per_cell);
 
  // check if sizes match
  AssertDimension(dof_values.size(), all_out.size());
 
  // check if we have enough dofs provided by the FE object
  // to interpolate the transferred data correctly
  for (auto it_dof_values = dof_values.begin();
       it_dof_values != dof_values.end();
       ++it_dof_values)
    Assert(
      dofs_per_cell == it_dof_values->size(),
      ExcMessage(
        "The transferred data was packed with a different number of dofs than the "
        "currently registered FE object assigned to the DoFHandler has."));
 
  // distribute data for each registered vector on mesh
  auto it_input  = dof_values.cbegin();
  auto it_output = all_out.begin();
  for (; it_input != dof_values.cend(); ++it_input, ++it_output)
    if (average_values)
      cell->distribute_local_to_global_by_interpolation(*it_input,
                                                        *(*it_output),
                                                        fe_index);
    else
      cell->set_dof_values_by_interpolation(*it_input,
                                            *(*it_output),
                                            fe_index,
                                            true);
 
  if (average_values)
    {
      // compute valence vector if averaging should be performed
      Vector<typename VectorType::value_type> ones(dofs_per_cell);
      ones = 1.0;
      cell->distribute_local_to_global_by_interpolation(ones,
                                                        valence,
                                                        fe_index);
    }
}
 
 
 
template <int dim, typename VectorType, int spacedim>
void
SolutionTransfer<dim, VectorType, spacedim>::clear()
{
  // nothing to do
}
 
 
 
namespace Legacy
{
 
  template <int dim, typename VectorType, int spacedim>
  SolutionTransfer<dim, VectorType, spacedim>::SolutionTransfer(
    const DoFHandler<dim, spacedim> &dof)
    : dof_handler(&dof, typeid(*this).name())
    , n_dofs_old(0)
    , prepared_for(none)
  {
    Assert((dynamic_cast<
              const parallel::distributed::Triangulation<dim, spacedim> *>(
              &dof_handler->get_triangulation()) == nullptr),
           ExcMessage(
             "You are calling the ::SolutionTransfer class "
             "with a DoFHandler that is built on a "
             "parallel::distributed::Triangulation. This will not "
             "work for parallel computations. You probably want to "
             "use the parallel::distributed::SolutionTransfer class."));
  }
 
 
 
  template <int dim, typename VectorType, int spacedim>
  SolutionTransfer<dim, VectorType, spacedim>::~SolutionTransfer()
  {
    clear();
  }
 
 
 
  template <int dim, typename VectorType, int spacedim>
  void
  SolutionTransfer<dim, VectorType, spacedim>::clear()
  {
    indices_on_cell.clear();
    dof_values_on_cell.clear();
    cell_map.clear();
 
    prepared_for = none;
  }
 
 
 
  template <int dim, typename VectorType, int spacedim>
  void
  SolutionTransfer<dim, VectorType, spacedim>::prepare_for_pure_refinement()
  {
    Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
    Assert(prepared_for != coarsening_and_refinement,
           ExcAlreadyPrepForCoarseAndRef());
 
    clear();
 
    // We need to access dof indices on the entire domain. For
    // parallel::shared::Triangulations, ownership of cells might change. If
    // they allow artificial cells, we need to restore the "true" cell owners
    // temporarily.
    // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
    // the current set of subdomain ids, set subdomain ids to the "true" owner
    // of each cell upon construction of the TemporarilyRestoreSubdomainIds
    // object, and later restore these flags when it is destroyed.
    const ::internal::parallel::shared::
      TemporarilyRestoreSubdomainIds<dim, spacedim>
        subdomain_modifier(dof_handler->get_triangulation());
 
    const unsigned int n_active_cells =
      dof_handler->get_triangulation().n_active_cells();
    n_dofs_old = dof_handler->n_dofs();
 
    // efficient reallocation of indices_on_cell
    std::vector<std::vector<types::global_dof_index>>(n_active_cells)
      .swap(indices_on_cell);
 
    for (const auto &cell : dof_handler->active_cell_iterators())
      {
        const unsigned int i = cell->active_cell_index();
        indices_on_cell[i].resize(cell->get_fe().n_dofs_per_cell());
        // on each cell store the indices of the
        // dofs. after refining we get the values
        // on the children by taking these
        // indices, getting the respective values
        // out of the data vectors and prolonging
        // them to the children
        cell->get_dof_indices(indices_on_cell[i]);
        cell_map[std::make_pair(cell->level(), cell->index())] =
          Pointerstruct(&indices_on_cell[i], cell->active_fe_index());
      }
    prepared_for = pure_refinement;
  }
 
 
 
  template <int dim, typename VectorType, int spacedim>
  void
  SolutionTransfer<dim, VectorType, spacedim>::refine_interpolate(
    const VectorType &in,
    VectorType       &out) const
  {
    Assert(prepared_for == pure_refinement, ExcNotPrepared());
    Assert(in.size() == n_dofs_old,
           ExcDimensionMismatch(in.size(), n_dofs_old));
    Assert(out.size() == dof_handler->n_dofs(),
           ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
    Assert(&in != &out,
           ExcMessage("Vectors cannot be used as input and output"
                      " at the same time!"));
 
    // We need to access dof indices on the entire domain. For
    // parallel::shared::Triangulations, ownership of cells might change. If
    // they allow artificial cells, we need to restore the "true" cell owners
    // temporarily.
    // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
    // the current set of subdomain ids, set subdomain ids to the "true" owner
    // of each cell upon construction of the TemporarilyRestoreSubdomainIds
    // object, and later restore these flags when it is destroyed.
    const ::internal::parallel::shared::
      TemporarilyRestoreSubdomainIds<dim, spacedim>
        subdomain_modifier(dof_handler->get_triangulation());
 
    Vector<typename VectorType::value_type> local_values(0);
 
    typename std::map<std::pair<unsigned int, unsigned int>,
                      Pointerstruct>::const_iterator pointerstruct,
      cell_map_end = cell_map.end();
 
    for (const auto &cell : dof_handler->cell_iterators())
      {
        pointerstruct =
          cell_map.find(std::make_pair(cell->level(), cell->index()));
 
        if (pointerstruct != cell_map_end)
          // this cell was refined or not
          // touched at all, so we can get
          // the new values by just setting
          // or interpolating to the children,
          // which is both done by one
          // function
          {
            const unsigned int this_fe_index =
              pointerstruct->second.active_fe_index;
            const unsigned int dofs_per_cell =
              cell->get_dof_handler().get_fe(this_fe_index).n_dofs_per_cell();
            local_values.reinit(dofs_per_cell, true);
 
            // make sure that the size of the stored indices is the same as
            // dofs_per_cell. since we store the desired fe_index, we know
            // what this size should be
            Assert(dofs_per_cell == (*pointerstruct->second.indices_ptr).size(),
                   ExcInternalError());
            for (unsigned int i = 0; i < dofs_per_cell; ++i)
              local_values(i) =
                ::internal::ElementAccess<VectorType>::get(
                  in, (*pointerstruct->second.indices_ptr)[i]);
            cell->set_dof_values_by_interpolation(local_values,
                                                  out,
                                                  this_fe_index,
                                                  true);
          }
      }
  }
 
 
 
  namespace internal
  {
    template <int dim, int spacedim>
    void
    extract_interpolation_matrices(
      const DoFHandler<dim, spacedim>      &dof,
      ::Table<2, FullMatrix<double>> &matrices)
    {
      if (dof.has_hp_capabilities() == false)
        return;
 
      const ::hp::FECollection<dim, spacedim> &fe =
        dof.get_fe_collection();
      matrices.reinit(fe.size(), fe.size());
      for (unsigned int i = 0; i < fe.size(); ++i)
        for (unsigned int j = 0; j < fe.size(); ++j)
          if (i != j)
            {
              matrices(i, j).reinit(fe[i].n_dofs_per_cell(),
                                    fe[j].n_dofs_per_cell());
 
              // see if we can get the interpolation matrices for this
              // combination of elements. if not, reset the matrix sizes to zero
              // to indicate that this particular combination isn't
              // supported. this isn't an outright error right away since we may
              // never need to actually interpolate between these two elements
              // on actual cells; we simply have to trigger an error if someone
              // actually tries
              try
                {
                  fe[i].get_interpolation_matrix(fe[j], matrices(i, j));
                }
              catch (const typename FiniteElement<dim, spacedim>::
                       ExcInterpolationNotImplemented &)
                {
                  matrices(i, j).reinit(0, 0);
                }
            }
    }
 
 
    template <int dim, int spacedim>
    void
    restriction_additive(const FiniteElement<dim, spacedim> &,
                         std::vector<std::vector<bool>> &)
    {}
 
    template <int dim, int spacedim>
    void
    restriction_additive(
      const ::hp::FECollection<dim, spacedim> &fe,
      std::vector<std::vector<bool>>                &restriction_is_additive)
    {
      restriction_is_additive.resize(fe.size());
      for (unsigned int f = 0; f < fe.size(); ++f)
        {
          restriction_is_additive[f].resize(fe[f].n_dofs_per_cell());
          for (unsigned int i = 0; i < fe[f].n_dofs_per_cell(); ++i)
            restriction_is_additive[f][i] = fe[f].restriction_is_additive(i);
        }
    }
  } // namespace internal
 
 
 
  template <int dim, typename VectorType, int spacedim>
  void
  SolutionTransfer<dim, VectorType, spacedim>::
    prepare_for_coarsening_and_refinement(const std::vector<VectorType> &all_in)
  {
    Assert(prepared_for != pure_refinement, ExcAlreadyPrepForRef());
    Assert(prepared_for != coarsening_and_refinement,
           ExcAlreadyPrepForCoarseAndRef());
 
    clear();
    n_dofs_old                 = dof_handler->n_dofs();
    const unsigned int in_size = all_in.size();
 
#ifdef DEBUG
    Assert(in_size != 0,
           ExcMessage("The array of input vectors you pass to this "
                      "function has no elements. This is not useful."));
    for (unsigned int i = 0; i < in_size; ++i)
      {
        Assert(all_in[i].size() == n_dofs_old,
               ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
      }
#endif
 
    // We need to access dof indices on the entire domain. For
    // parallel::shared::Triangulations, ownership of cells might change. If
    // they allow artificial cells, we need to restore the "true" cell owners
    // temporarily.
    // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
    // the current set of subdomain ids, set subdomain ids to the "true" owner
    // of each cell upon construction of the TemporarilyRestoreSubdomainIds
    // object, and later restore these flags when it is destroyed.
    const ::internal::parallel::shared::
      TemporarilyRestoreSubdomainIds<dim, spacedim>
        subdomain_modifier(dof_handler->get_triangulation());
 
    // first count the number
    // of cells that will be coarsened
    // and that'll stay or be refined
    unsigned int n_cells_to_coarsen        = 0;
    unsigned int n_cells_to_stay_or_refine = 0;
    for (const auto &act_cell : dof_handler->active_cell_iterators())
      {
        if (act_cell->coarsen_flag_set())
          ++n_cells_to_coarsen;
        else
          ++n_cells_to_stay_or_refine;
      }
    Assert((n_cells_to_coarsen + n_cells_to_stay_or_refine) ==
             dof_handler->get_triangulation().n_active_cells(),
           ExcInternalError());
 
    unsigned int n_coarsen_fathers = 0;
    for (const auto &cell : dof_handler->cell_iterators())
      if (!cell->is_active() && cell->child(0)->coarsen_flag_set())
        ++n_coarsen_fathers;
    Assert(n_cells_to_coarsen >= 2 * n_coarsen_fathers, ExcInternalError());
    (void)n_cells_to_coarsen;
 
    // allocate the needed memory. initialize
    // the following arrays in an efficient
    // way, without copying much
    std::vector<std::vector<types::global_dof_index>>(n_cells_to_stay_or_refine)
      .swap(indices_on_cell);
 
    std::vector<std::vector<Vector<typename VectorType::value_type>>>(
      n_coarsen_fathers,
      std::vector<Vector<typename VectorType::value_type>>(in_size))
      .swap(dof_values_on_cell);
 
    Table<2, FullMatrix<double>>   interpolation_hp;
    std::vector<std::vector<bool>> restriction_is_additive;
 
    internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
    internal::restriction_additive(dof_handler->get_fe_collection(),
                                   restriction_is_additive);
 
    // we need counters for
    // the 'to_stay_or_refine' cells 'n_sr' and
    // the 'coarsen_fathers' cells 'n_cf',
    unsigned int n_sr = 0, n_cf = 0;
    for (const auto &cell : dof_handler->cell_iterators())
      {
        // CASE 1: active cell that remains as it is
        if (cell->is_active() && !cell->coarsen_flag_set())
          {
            const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
            indices_on_cell[n_sr].resize(dofs_per_cell);
            // cell will not be coarsened,
            // so we get away by storing the
            // dof indices and later
            // interpolating to the children
            cell->get_dof_indices(indices_on_cell[n_sr]);
            cell_map[std::make_pair(cell->level(), cell->index())] =
              Pointerstruct(&indices_on_cell[n_sr], cell->active_fe_index());
            ++n_sr;
          }
 
        // CASE 2: cell is inactive but will become active
        else if (cell->has_children() && cell->child(0)->coarsen_flag_set())
          {
            // we will need to interpolate from the children of this cell
            // to the current one. in the hp-context, this also means
            // we need to figure out which finite element space to interpolate
            // to since that is not implied by the global FE as in the non-hp-
            // case. we choose the 'least dominant fe' on all children from
            // the associated FECollection.
            std::set<unsigned int> fe_indices_children;
            for (const auto &child : cell->child_iterators())
              {
                Assert(child->is_active() && child->coarsen_flag_set(),
                       typename ::Triangulation<
                         dim>::ExcInconsistentCoarseningFlags());
 
                fe_indices_children.insert(child->active_fe_index());
              }
            Assert(!fe_indices_children.empty(), ExcInternalError());
 
            const unsigned int target_fe_index =
              dof_handler->get_fe_collection().find_dominated_fe_extended(
                fe_indices_children, /*codim=*/0);
 
            Assert(target_fe_index != numbers::invalid_unsigned_int,
                   ::internal::hp::DoFHandlerImplementation::
                     ExcNoDominatedFiniteElementOnChildren());
 
            const unsigned int dofs_per_cell =
              dof_handler->get_fe(target_fe_index).n_dofs_per_cell();
 
            std::vector<Vector<typename VectorType::value_type>>(
              in_size, Vector<typename VectorType::value_type>(dofs_per_cell))
              .swap(dof_values_on_cell[n_cf]);
 
 
            // store the data of each of the input vectors. get this data
            // as interpolated onto a finite element space that encompasses
            // that of all the children. note that
            // cell->get_interpolated_dof_values already does all of the
            // interpolations between spaces
            for (unsigned int j = 0; j < in_size; ++j)
              cell->get_interpolated_dof_values(all_in[j],
                                                dof_values_on_cell[n_cf][j],
                                                target_fe_index);
            cell_map[std::make_pair(cell->level(), cell->index())] =
              Pointerstruct(&dof_values_on_cell[n_cf], target_fe_index);
            ++n_cf;
          }
      }
    Assert(n_sr == n_cells_to_stay_or_refine, ExcInternalError());
    Assert(n_cf == n_coarsen_fathers, ExcInternalError());
 
    prepared_for = coarsening_and_refinement;
  }
 
 
 
  template <int dim, typename VectorType, int spacedim>
  void
  SolutionTransfer<dim, VectorType, spacedim>::
    prepare_for_coarsening_and_refinement(const VectorType &in)
  {
    std::vector<VectorType> all_in(1, in);
    prepare_for_coarsening_and_refinement(all_in);
  }
 
 
 
  template <int dim, typename VectorType, int spacedim>
  void
  SolutionTransfer<dim, VectorType, spacedim>::interpolate(
    const std::vector<VectorType> &all_in,
    std::vector<VectorType>       &all_out) const
  {
    const unsigned int size = all_in.size();
#ifdef DEBUG
    Assert(prepared_for == coarsening_and_refinement, ExcNotPrepared());
    Assert(all_out.size() == size, ExcDimensionMismatch(all_out.size(), size));
    for (unsigned int i = 0; i < size; ++i)
      Assert(all_in[i].size() == n_dofs_old,
             ExcDimensionMismatch(all_in[i].size(), n_dofs_old));
    for (unsigned int i = 0; i < all_out.size(); ++i)
      Assert(all_out[i].size() == dof_handler->n_dofs(),
             ExcDimensionMismatch(all_out[i].size(), dof_handler->n_dofs()));
    for (unsigned int i = 0; i < size; ++i)
      for (unsigned int j = 0; j < size; ++j)
        Assert(&all_in[i] != &all_out[j],
               ExcMessage("Vectors cannot be used as input and output"
                          " at the same time!"));
#endif
 
    // We need to access dof indices on the entire domain. For
    // parallel::shared::Triangulations, ownership of cells might change. If
    // they allow artificial cells, we need to restore the "true" cell owners
    // temporarily.
    // We use the TemporarilyRestoreSubdomainIds class for this purpose: we save
    // the current set of subdomain ids, set subdomain ids to the "true" owner
    // of each cell upon construction of the TemporarilyRestoreSubdomainIds
    // object, and later restore these flags when it is destroyed.
    const ::internal::parallel::shared::
      TemporarilyRestoreSubdomainIds<dim, spacedim>
        subdomain_modifier(dof_handler->get_triangulation());
 
    Vector<typename VectorType::value_type> local_values;
    std::vector<types::global_dof_index>    dofs;
 
    typename std::map<std::pair<unsigned int, unsigned int>,
                      Pointerstruct>::const_iterator pointerstruct,
      cell_map_end = cell_map.end();
 
    Table<2, FullMatrix<double>> interpolation_hp;
    internal::extract_interpolation_matrices(*dof_handler, interpolation_hp);
    Vector<typename VectorType::value_type> tmp, tmp2;
 
    for (const auto &cell : dof_handler->cell_iterators())
      {
        pointerstruct =
          cell_map.find(std::make_pair(cell->level(), cell->index()));
 
        if (pointerstruct != cell_map_end)
          {
            const std::vector<types::global_dof_index> *const indexptr =
              pointerstruct->second.indices_ptr;
 
            const std::vector<Vector<typename VectorType::value_type>>
              *const valuesptr = pointerstruct->second.dof_values_ptr;
 
            // cell stayed as it was or was refined
            if (indexptr != nullptr)
              {
                Assert(valuesptr == nullptr, ExcInternalError());
 
                const unsigned int old_fe_index =
                  pointerstruct->second.active_fe_index;
 
                // get the values of each of the input data vectors on this cell
                // and prolong it to its children
                unsigned int in_size = indexptr->size();
                for (unsigned int j = 0; j < size; ++j)
                  {
                    tmp.reinit(in_size, true);
                    for (unsigned int i = 0; i < in_size; ++i)
                      tmp(i) = ::internal::ElementAccess<VectorType>::get(
                        all_in[j], (*indexptr)[i]);
 
                    cell->set_dof_values_by_interpolation(tmp,
                                                          all_out[j],
                                                          old_fe_index,
                                                          true);
                  }
              }
            else if (valuesptr)
              // the children of this cell were deleted
              {
                Assert(!cell->has_children(), ExcInternalError());
                Assert(indexptr == nullptr, ExcInternalError());
 
                const unsigned int dofs_per_cell =
                  cell->get_fe().n_dofs_per_cell();
                dofs.resize(dofs_per_cell);
                // get the local
                // indices
                cell->get_dof_indices(dofs);
 
                // distribute the stored data to the new vectors
                for (unsigned int j = 0; j < size; ++j)
                  {
                    // make sure that the size of the stored indices is the same
                    // as dofs_per_cell. this is kind of a test if we use the
                    // same FE in the hp-case. to really do that test we would
                    // have to store the fe_index of all cells
                    const Vector<typename VectorType::value_type> *data =
                      nullptr;
                    const unsigned int active_fe_index =
                      cell->active_fe_index();
                    if (active_fe_index !=
                        pointerstruct->second.active_fe_index)
                      {
                        const unsigned int old_index =
                          pointerstruct->second.active_fe_index;
                        const FullMatrix<double> &interpolation_matrix =
                          interpolation_hp(active_fe_index, old_index);
                        // The interpolation matrix might be empty when using
                        // FE_Nothing.
                        if (interpolation_matrix.empty())
                          tmp.reinit(dofs_per_cell, false);
                        else
                          {
                            tmp.reinit(dofs_per_cell, true);
                            AssertDimension((*valuesptr)[j].size(),
                                            interpolation_matrix.n());
                            AssertDimension(tmp.size(),
                                            interpolation_matrix.m());
                            interpolation_matrix.vmult(tmp, (*valuesptr)[j]);
                          }
                        data = &tmp;
                      }
                    else
                      data = &(*valuesptr)[j];
 
 
                    for (unsigned int i = 0; i < dofs_per_cell; ++i)
                      ::internal::ElementAccess<VectorType>::set(
                        (*data)(i), dofs[i], all_out[j]);
                  }
              }
            // undefined status
            else
              Assert(false, ExcInternalError());
          }
      }
 
    // We have written into the output vectors. If this was a PETSc vector, for
    // example, then we need to compress these to make future operations safe:
    for (auto &vec : all_out)
      vec.compress(VectorOperation::insert);
  }
 
 
 
  template <int dim, typename VectorType, int spacedim>
  void
  SolutionTransfer<dim, VectorType, spacedim>::interpolate(
    const VectorType &in,
    VectorType       &out) const
  {
    Assert(in.size() == n_dofs_old,
           ExcDimensionMismatch(in.size(), n_dofs_old));
    Assert(out.size() == dof_handler->n_dofs(),
           ExcDimensionMismatch(out.size(), dof_handler->n_dofs()));
 
    std::vector<VectorType> all_in  = {in};
    std::vector<VectorType> all_out = {out};
 
    interpolate(all_in, all_out);
 
    out = all_out[0];
  }
 
 
 
  template <int dim, typename VectorType, int spacedim>
  std::size_t
  SolutionTransfer<dim, VectorType, spacedim>::memory_consumption() const
  {
    // at the moment we do not include the memory
    // consumption of the cell_map as we have no
    // real idea about memory consumption of a
    // std::map
    return (MemoryConsumption::memory_consumption(dof_handler) +
            MemoryConsumption::memory_consumption(n_dofs_old) +
            sizeof(prepared_for) +
            MemoryConsumption::memory_consumption(indices_on_cell) +
            MemoryConsumption::memory_consumption(dof_values_on_cell));
  }
 
 
 
  template <int dim, typename VectorType, int spacedim>
  std::size_t
  SolutionTransfer<dim, VectorType, spacedim>::Pointerstruct::
    memory_consumption() const
  {
    return sizeof(*this);
  }
 
} // namespace Legacy
 
 
/*-------------- Explicit Instantiations -------------------------------*/
#define SPLIT_INSTANTIATIONS_COUNT 4
#ifndef SPLIT_INSTANTIATIONS_INDEX
#  define SPLIT_INSTANTIATIONS_INDEX 0
#endif
#include "solution_transfer.inst"
 
DEAL_II_NAMESPACE_CLOSE