dof_handler_policy.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


#include <deal.II/base/geometry_info.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/partitioner.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/work_stream.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_handler_policy.h>

#include <deal.II/fe/fe.h>

#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator.h>

#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#ifdef DEAL_II_WITH_ZLIB
#  include <boost/iostreams/device/back_inserter.hpp>
#  include <boost/iostreams/filter/gzip.hpp>
#  include <boost/iostreams/filtering_stream.hpp>
#  include <boost/iostreams/stream.hpp>
#  include <boost/serialization/array.hpp>
#endif

#include <algorithm>
#include <memory>
#include <numeric>
#include <set>

DEAL_II_NAMESPACE_OPEN


namespace internal
{
  namespace DoFHandlerImplementation
  {
    namespace Policy
    {
      // use class ::DoFHandler instead
      // of namespace internal::DoFHandler in
      // the following
      using ::DoFHandler;

      namespace hp
      {
        using ::hp::DoFHandler;
      }


      namespace
      {
        template <class DoFHandlerType>
        void
        update_all_active_cell_dof_indices_caches(
          const DoFHandlerType &dof_handler)
        {
          typename DoFHandlerType::active_cell_iterator
            beginc = dof_handler.begin_active(),
            endc   = dof_handler.end();

          auto worker =
            [](const typename DoFHandlerType::active_cell_iterator &cell,
               void *,
               void *) {
              if (!cell->is_artificial())
                cell->update_cell_dof_indices_cache();
            };

          // parallelize filling all of the cell caches. by using
          // WorkStream, we make sure that we only run through the
          // range of iterators once, whereas a parallel_for loop
          // for example has to split the range multiple times,
          // which is expensive because cell iterators are not
          // random access iterators with a cheap operator-
          WorkStream::run(beginc,
                          endc,
                          worker,
                          /* copier */ std::function<void(void *)>(),
                          /* scratch_data */ nullptr,
                          /* copy_data */ nullptr,
                          2 * MultithreadInfo::n_threads(),
                          /* chunk_size = */ 32);
        }


        template <class DoFHandlerType>
        void
        update_all_level_cell_dof_indices_caches(
          const DoFHandlerType &dof_handler)
        {
          typename DoFHandlerType::level_cell_iterator beginc =
                                                         dof_handler.begin(),
                                                       endc = dof_handler.end();

          auto worker =
            [](const typename DoFHandlerType::level_cell_iterator &cell,
               void *,
               void *) {
              if (cell->has_children() || !cell->is_artificial())
                cell->update_cell_dof_indices_cache();
            };

          // parallelize filling all of the cell caches. by using
          // WorkStream, we make sure that we only run through the
          // range of iterators once, whereas a parallel_for loop
          // for example has to split the range multiple times,
          // which is expensive because cell iterators are not
          // random access iterators with a cheap operator-
          WorkStream::run(beginc,
                          endc,
                          worker,
                          /* copier */ std::function<void(void *)>(),
                          /* scratch_data */ nullptr,
                          /* copy_data */ nullptr,
                          2 * MultithreadInfo::n_threads(),
                          /* chunk_size = */ 32);
        }


        using DoFIdentities =
          std::vector<std::pair<unsigned int, unsigned int>>;


        template <int structdim, int dim, int spacedim>
        void
        ensure_existence_of_dof_identities(
          const FiniteElement<dim, spacedim> &fe1,
          const FiniteElement<dim, spacedim> &fe2,
          std::unique_ptr<DoFIdentities> &    identities)
        {
          // see if we need to fill this entry, or whether it already
          // exists
          if (identities.get() == nullptr)
            {
              switch (structdim)
                {
                  case 0:
                    {
                      identities = std_cxx14::make_unique<DoFIdentities>(
                        fe1.hp_vertex_dof_identities(fe2));
                      break;
                    }

                  case 1:
                    {
                      identities = std_cxx14::make_unique<DoFIdentities>(
                        fe1.hp_line_dof_identities(fe2));
                      break;
                    }

                  case 2:
                    {
                      identities = std_cxx14::make_unique<DoFIdentities>(
                        fe1.hp_quad_dof_identities(fe2));
                      break;
                    }

                  default:
                    Assert(false, ExcNotImplemented());
                }

              // double check whether the newly created entries make
              // any sense at all
              for (unsigned int i = 0; i < identities->size(); ++i)
                {
                  Assert((*identities)[i].first <
                           fe1.template n_dofs_per_object<structdim>(),
                         ExcInternalError());
                  Assert((*identities)[i].second <
                           fe2.template n_dofs_per_object<structdim>(),
                         ExcInternalError());
                }
            }
        }



        template <int dim, int spacedim, typename iterator>
        unsigned int
        get_most_dominating_fe_index(const iterator &object)
        {
          unsigned int dominating_fe_index = 0;
          for (; dominating_fe_index < object->n_active_fe_indices();
               ++dominating_fe_index)
            {
              const FiniteElement<dim, spacedim> &this_fe = object->get_fe(
                object->nth_active_fe_index(dominating_fe_index));

              FiniteElementDomination::Domination domination =
                FiniteElementDomination::either_element_can_dominate;
              for (unsigned int other_fe_index = 0;
                   other_fe_index < object->n_active_fe_indices();
                   ++other_fe_index)
                if (other_fe_index != dominating_fe_index)
                  {
                    const FiniteElement<dim, spacedim> &that_fe =
                      object->get_fe(
                        object->nth_active_fe_index(other_fe_index));

                    domination =
                      domination & this_fe.compare_for_face_domination(that_fe);
                  }

              // see if this element is able to dominate all the other
              // ones, and if so take it
              if ((domination ==
                   FiniteElementDomination::this_element_dominates) ||
                  (domination ==
                   FiniteElementDomination::either_element_can_dominate) ||
                  (domination == FiniteElementDomination::no_requirements))
                break;
            }

          // check that we have found one such fe
          if (dominating_fe_index != object->n_active_fe_indices())
            {
              // return the finite element index used on it. note that
              // only a single fe can be active on such subfaces
              return object->nth_active_fe_index(dominating_fe_index);
            }
          else
            {
              // if we couldn't find the most dominating object
              return numbers::invalid_unsigned_int;
            }
        }
      } // namespace



      struct Implementation
      {
        /* -------------- distribute_dofs functionality ------------- */

        template <int spacedim>
        static types::global_dof_index
        distribute_dofs_on_cell(
          const DoFHandler<1, spacedim> &dof_handler,
          const typename DoFHandler<1, spacedim>::active_cell_iterator &cell,
          types::global_dof_index next_free_dof)
        {
          // distribute dofs of vertices
          if (dof_handler.get_fe().dofs_per_vertex > 0)
            for (unsigned int v = 0; v < GeometryInfo<1>::vertices_per_cell;
                 ++v)
              {
                if (cell->vertex_dof_index(v, 0) == numbers::invalid_dof_index)
                  for (unsigned int d = 0;
                       d < dof_handler.get_fe().dofs_per_vertex;
                       ++d)
                    {
                      Assert((cell->vertex_dof_index(v, d) ==
                              numbers::invalid_dof_index),
                             ExcInternalError());
                      cell->set_vertex_dof_index(v, d, next_free_dof++);
                    }
                else
                  for (unsigned int d = 0;
                       d < dof_handler.get_fe().dofs_per_vertex;
                       ++d)
                    Assert((cell->vertex_dof_index(v, d) !=
                            numbers::invalid_dof_index),
                           ExcInternalError());
              }

          // dofs of line
          for (unsigned int d = 0; d < dof_handler.get_fe().dofs_per_line; ++d)
            cell->set_dof_index(d, next_free_dof++);

          return next_free_dof;
        }



        template <int spacedim>
        static types::global_dof_index
        distribute_dofs_on_cell(
          const DoFHandler<2, spacedim> &dof_handler,
          const typename DoFHandler<2, spacedim>::active_cell_iterator &cell,
          types::global_dof_index next_free_dof)
        {
          if (dof_handler.get_fe().dofs_per_vertex > 0)
            // number dofs on vertices
            for (unsigned int vertex = 0;
                 vertex < GeometryInfo<2>::vertices_per_cell;
                 ++vertex)
              // check whether dofs for this vertex have been distributed
              // (checking the first dof should be good enough)
              if (cell->vertex_dof_index(vertex, 0) ==
                  numbers::invalid_dof_index)
                for (unsigned int d = 0;
                     d < dof_handler.get_fe().dofs_per_vertex;
                     ++d)
                  cell->set_vertex_dof_index(vertex, d, next_free_dof++);

          // for the four sides
          if (dof_handler.get_fe().dofs_per_line > 0)
            for (unsigned int side = 0; side < GeometryInfo<2>::faces_per_cell;
                 ++side)
              {
                const typename DoFHandler<2, spacedim>::line_iterator line =
                  cell->line(side);

                // distribute dofs if necessary: check whether line dof is
                // already numbered (checking the first dof should be good
                // enough)
                if (line->dof_index(0) == numbers::invalid_dof_index)
                  // if not: distribute dofs
                  for (unsigned int d = 0;
                       d < dof_handler.get_fe().dofs_per_line;
                       ++d)
                    line->set_dof_index(d, next_free_dof++);
              }


          // dofs of quad
          if (dof_handler.get_fe().dofs_per_quad > 0)
            for (unsigned int d = 0; d < dof_handler.get_fe().dofs_per_quad;
                 ++d)
              cell->set_dof_index(d, next_free_dof++);

          return next_free_dof;
        }



        template <int spacedim>
        static types::global_dof_index
        distribute_dofs_on_cell(
          const DoFHandler<3, spacedim> &dof_handler,
          const typename DoFHandler<3, spacedim>::active_cell_iterator &cell,
          types::global_dof_index next_free_dof)
        {
          if (dof_handler.get_fe().dofs_per_vertex > 0)
            // number dofs on vertices
            for (unsigned int vertex = 0;
                 vertex < GeometryInfo<3>::vertices_per_cell;
                 ++vertex)
              // check whether dofs for this vertex have been distributed
              // (checking the first dof should be good enough)
              if (cell->vertex_dof_index(vertex, 0) ==
                  numbers::invalid_dof_index)
                for (unsigned int d = 0;
                     d < dof_handler.get_fe().dofs_per_vertex;
                     ++d)
                  cell->set_vertex_dof_index(vertex, d, next_free_dof++);

          // for the lines
          if (dof_handler.get_fe().dofs_per_line > 0)
            for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell; ++l)
              {
                const typename DoFHandler<3, spacedim>::line_iterator line =
                  cell->line(l);

                // distribute dofs if necessary: check whether line dof is
                // already numbered (checking the first dof should be good
                // enough)
                if (line->dof_index(0) == numbers::invalid_dof_index)
                  // if not: distribute dofs
                  for (unsigned int d = 0;
                       d < dof_handler.get_fe().dofs_per_line;
                       ++d)
                    line->set_dof_index(d, next_free_dof++);
              }

          // for the quads
          if (dof_handler.get_fe().dofs_per_quad > 0)
            for (unsigned int q = 0; q < GeometryInfo<3>::quads_per_cell; ++q)
              {
                const typename DoFHandler<3, spacedim>::quad_iterator quad =
                  cell->quad(q);

                // distribute dofs if necessary: check whether line dof is
                // already numbered (checking the first dof should be good
                // enough)
                if (quad->dof_index(0) == numbers::invalid_dof_index)
                  // if not: distribute dofs
                  for (unsigned int d = 0;
                       d < dof_handler.get_fe().dofs_per_quad;
                       ++d)
                    quad->set_dof_index(d, next_free_dof++);
              }


          // dofs of hex
          if (dof_handler.get_fe().dofs_per_hex > 0)
            for (unsigned int d = 0; d < dof_handler.get_fe().dofs_per_hex; ++d)
              cell->set_dof_index(d, next_free_dof++);

          return next_free_dof;
        }



        // same for the hp::DoFHandler
        template <int spacedim>
        static types::global_dof_index
        distribute_dofs_on_cell(
          const hp::DoFHandler<1, spacedim> &,
          const typename hp::DoFHandler<1, spacedim>::active_cell_iterator
            &                     cell,
          types::global_dof_index next_free_dof)
        {
          const unsigned int dim = 1;

          const FiniteElement<dim, spacedim> &fe = cell->get_fe();
          const unsigned int fe_index            = cell->active_fe_index();

          // number dofs on vertices. to do so, check whether dofs for
          // this vertex have been distributed and for the present fe
          // (only check the first dof), and if this isn't the case
          // distribute new ones there
          if (fe.dofs_per_vertex > 0)
            for (unsigned int vertex = 0;
                 vertex < GeometryInfo<1>::vertices_per_cell;
                 ++vertex)
              if (cell->vertex_dof_index(vertex, 0, fe_index) ==
                  numbers::invalid_dof_index)
                for (unsigned int d = 0; d < fe.dofs_per_vertex;
                     ++d, ++next_free_dof)
                  cell->set_vertex_dof_index(vertex,
                                             d,
                                             next_free_dof,
                                             fe_index);

          // finally for the line. this one shouldn't be numbered yet
          if (fe.dofs_per_line > 0)
            {
              Assert((cell->dof_index(0, fe_index) ==
                      numbers::invalid_dof_index),
                     ExcInternalError());

              for (unsigned int d = 0; d < fe.dofs_per_line;
                   ++d, ++next_free_dof)
                cell->set_dof_index(d, next_free_dof, fe_index);
            }

          // note that this cell has been processed
          cell->set_user_flag();

          return next_free_dof;
        }



        template <int spacedim>
        static types::global_dof_index
        distribute_dofs_on_cell(
          const hp::DoFHandler<2, spacedim> &,
          const typename hp::DoFHandler<2, spacedim>::active_cell_iterator
            &                     cell,
          types::global_dof_index next_free_dof)
        {
          const unsigned int dim = 2;

          const FiniteElement<dim, spacedim> &fe = cell->get_fe();
          const unsigned int fe_index            = cell->active_fe_index();

          // number dofs on vertices. to do so, check whether dofs for
          // this vertex have been distributed and for the present fe
          // (only check the first dof), and if this isn't the case
          // distribute new ones there
          if (fe.dofs_per_vertex > 0)
            for (unsigned int vertex = 0;
                 vertex < GeometryInfo<2>::vertices_per_cell;
                 ++vertex)
              if (cell->vertex_dof_index(vertex, 0, fe_index) ==
                  numbers::invalid_dof_index)
                for (unsigned int d = 0; d < fe.dofs_per_vertex;
                     ++d, ++next_free_dof)
                  cell->set_vertex_dof_index(vertex,
                                             d,
                                             next_free_dof,
                                             fe_index);

          // next the sides. do the same as above: check whether the
          // line is already numbered for the present fe_index, and if
          // not do it
          if (fe.dofs_per_line > 0)
            for (unsigned int l = 0; l < GeometryInfo<2>::lines_per_cell; ++l)
              {
                typename hp::DoFHandler<dim, spacedim>::line_iterator line =
                  cell->line(l);

                if (line->dof_index(0, fe_index) == numbers::invalid_dof_index)
                  for (unsigned int d = 0; d < fe.dofs_per_line;
                       ++d, ++next_free_dof)
                    line->set_dof_index(d, next_free_dof, fe_index);
              }


          // finally for the quad. this one shouldn't be numbered yet
          if (fe.dofs_per_quad > 0)
            {
              Assert((cell->dof_index(0, fe_index) ==
                      numbers::invalid_dof_index),
                     ExcInternalError());

              for (unsigned int d = 0; d < fe.dofs_per_quad;
                   ++d, ++next_free_dof)
                cell->set_dof_index(d, next_free_dof, fe_index);
            }

          // note that this cell has been processed
          cell->set_user_flag();

          return next_free_dof;
        }



        template <int spacedim>
        static types::global_dof_index
        distribute_dofs_on_cell(
          const hp::DoFHandler<3, spacedim> &,
          const typename hp::DoFHandler<3, spacedim>::active_cell_iterator
            &                     cell,
          types::global_dof_index next_free_dof)
        {
          const unsigned int dim = 3;

          const FiniteElement<dim, spacedim> &fe = cell->get_fe();
          const unsigned int fe_index            = cell->active_fe_index();

          // number dofs on vertices. to do so, check whether dofs for
          // this vertex have been distributed and for the present fe
          // (only check the first dof), and if this isn't the case
          // distribute new ones there
          if (fe.dofs_per_vertex > 0)
            for (unsigned int vertex = 0;
                 vertex < GeometryInfo<3>::vertices_per_cell;
                 ++vertex)
              if (cell->vertex_dof_index(vertex, 0, fe_index) ==
                  numbers::invalid_dof_index)
                for (unsigned int d = 0; d < fe.dofs_per_vertex;
                     ++d, ++next_free_dof)
                  cell->set_vertex_dof_index(vertex,
                                             d,
                                             next_free_dof,
                                             fe_index);

          // next the four lines. do the same as above: check whether
          // the line is already numbered for the present fe_index,
          // and if not do it
          if (fe.dofs_per_line > 0)
            for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell; ++l)
              {
                typename hp::DoFHandler<dim, spacedim>::line_iterator line =
                  cell->line(l);

                if (line->dof_index(0, fe_index) == numbers::invalid_dof_index)
                  for (unsigned int d = 0; d < fe.dofs_per_line;
                       ++d, ++next_free_dof)
                    line->set_dof_index(d, next_free_dof, fe_index);
              }

          // same for quads
          if (fe.dofs_per_quad > 0)
            for (unsigned int q = 0; q < GeometryInfo<3>::quads_per_cell; ++q)
              {
                typename hp::DoFHandler<dim, spacedim>::quad_iterator quad =
                  cell->quad(q);

                if (quad->dof_index(0, fe_index) == numbers::invalid_dof_index)
                  for (unsigned int d = 0; d < fe.dofs_per_quad;
                       ++d, ++next_free_dof)
                    quad->set_dof_index(d, next_free_dof, fe_index);
              }


          // finally for the hex. this one shouldn't be numbered yet
          if (fe.dofs_per_hex > 0)
            {
              Assert((cell->dof_index(0, fe_index) ==
                      numbers::invalid_dof_index),
                     ExcInternalError());

              for (unsigned int d = 0; d < fe.dofs_per_hex;
                   ++d, ++next_free_dof)
                cell->set_dof_index(d, next_free_dof, fe_index);
            }

          // note that this cell has been processed
          cell->set_user_flag();

          return next_free_dof;
        }



        template <int dim, int spacedim>
        static std::map<types::global_dof_index, types::global_dof_index>
        compute_vertex_dof_identities(
          hp::DoFHandler<dim, spacedim> &dof_handler)
        {
          std::map<types::global_dof_index, types::global_dof_index>
            dof_identities;

          // Note: we may wish to have something here similar to what
          // we do for lines and quads, namely that we only identify
          // dofs for any fe towards the most dominating one. however,
          // it is not clear whether this is actually necessary for
          // vertices at all, I can't think of a finite element that
          // would make that necessary...
          ::Table<2, std::unique_ptr<DoFIdentities>>
            vertex_dof_identities(dof_handler.get_fe_collection().size(),
                                  dof_handler.get_fe_collection().size());

          // first identify vertices we want to exclude from working on.
          // specifically, these are the vertices of artificial and ghost
          // cells because at the time when we get here, we do not yet
          // know DoF indices on ghost cells (and we will never know
          // them for artificial cells). this is, at least the case for
          // parallel::distributed::Triangulations.
          //
          // this means that we will not unify DoF indices between locally
          // owned cells and ghost cells, and this is different from what
          // we would do if the triangulation were not split into subdomains.
          // on the other hand, DoF unification is only an optimization: we
          // will still record these identities when we compute hanging
          // node constraints; we just end up with more DoFs than we would
          // if we unified DoF indices also between locally owned and ghost
          // cells, but we end up with a simpler algorithm in return.
          std::vector<bool> include_vertex =
            dof_handler.get_triangulation().get_used_vertices();
          if (dynamic_cast<
                const parallel::distributed::Triangulation<dim, spacedim> *>(
                &dof_handler.get_triangulation()) != nullptr)
            for (const auto &cell : dof_handler.active_cell_iterators())
              if (!cell->is_locally_owned())
                for (unsigned int v = 0;
                     v < GeometryInfo<dim>::vertices_per_cell;
                     ++v)
                  include_vertex[cell->vertex_index(v)] = false;

          // loop over all vertices and see which one we need to work
          // on
          for (unsigned int vertex_index = 0;
               vertex_index < dof_handler.get_triangulation().n_vertices();
               ++vertex_index)
            if ((dof_handler.get_triangulation()
                   .get_used_vertices()[vertex_index] == true) &&
                (include_vertex[vertex_index] == true))
              {
                const unsigned int n_active_fe_indices =
                  ::internal::DoFAccessorImplementation::Implementation::
                    n_active_vertex_fe_indices(dof_handler, vertex_index);
                if (n_active_fe_indices > 1)
                  {
                    const unsigned int first_fe_index =
                      ::internal::DoFAccessorImplementation::
                        Implementation::nth_active_vertex_fe_index(dof_handler,
                                                                   vertex_index,
                                                                   0);

                    // loop over all the other FEs with which we want
                    // to identify the DoF indices of the first FE of
                    for (unsigned int f = 1; f < n_active_fe_indices; ++f)
                      {
                        const unsigned int other_fe_index =
                          ::internal::DoFAccessorImplementation::
                            Implementation::nth_active_vertex_fe_index(
                              dof_handler, vertex_index, f);

                        // make sure the entry in the equivalence
                        // table exists
                        ensure_existence_of_dof_identities<0>(
                          dof_handler.get_fe(first_fe_index),
                          dof_handler.get_fe(other_fe_index),
                          vertex_dof_identities[first_fe_index]
                                               [other_fe_index]);

                        // then loop through the identities we
                        // have. first get the global numbers of the
                        // dofs we want to identify and make sure they
                        // are not yet constrained to anything else,
                        // except for to each other. use the rule that
                        // we will always constrain the dof with the
                        // higher fe index to the one with the lower,
                        // to avoid circular reasoning.
                        DoFIdentities &identities =
                          *vertex_dof_identities[first_fe_index]
                                                [other_fe_index];
                        for (unsigned int i = 0; i < identities.size(); ++i)
                          {
                            const types::global_dof_index lower_dof_index =
                              ::internal::DoFAccessorImplementation::
                                Implementation::get_vertex_dof_index(
                                  dof_handler,
                                  vertex_index,
                                  first_fe_index,
                                  identities[i].first);
                            const types::global_dof_index higher_dof_index =
                              ::internal::DoFAccessorImplementation::
                                Implementation::get_vertex_dof_index(
                                  dof_handler,
                                  vertex_index,
                                  other_fe_index,
                                  identities[i].second);

                            Assert((dof_identities.find(higher_dof_index) ==
                                    dof_identities.end()) ||
                                     (dof_identities[higher_dof_index] ==
                                      lower_dof_index),
                                   ExcInternalError());

                            dof_identities[higher_dof_index] = lower_dof_index;
                          }
                      }
                  }
              }

          return dof_identities;
        }


        template <int spacedim>
        static std::map<types::global_dof_index, types::global_dof_index>
          compute_line_dof_identities(hp::DoFHandler<1, spacedim> &)
        {
          return std::map<types::global_dof_index, types::global_dof_index>();
        }


        template <int dim, int spacedim>
        static std::map<types::global_dof_index, types::global_dof_index>
        compute_line_dof_identities(hp::DoFHandler<dim, spacedim> &dof_handler)
        {
          std::map<types::global_dof_index, types::global_dof_index>
            dof_identities;

          // we will mark lines that we have already treated, so first save and
          // clear the user flags on lines and later restore them
          std::vector<bool> user_flags;
          dof_handler.get_triangulation().save_user_flags_line(user_flags);
          const_cast<::Triangulation<dim, spacedim> &>(
            dof_handler.get_triangulation())
            .clear_user_flags_line();

          // exclude lines that bound cells we don't locally own, because
          // we do not have information about their dofs at this point.
          // this is, at least the case for
          // parallel::distributed::Triangulations.
          //
          // this means that we will not unify DoF indices between locally
          // owned cells and ghost cells, and this is different from what
          // we would do if the triangulation were not split into subdomains.
          // on the other hand, DoF unification is only an optimization: we
          // will still record these identities when we compute hanging
          // node constraints; we just end up with more DoFs than we would
          // if we unified DoF indices also between locally owned and ghost
          // cells, but we end up with a simpler algorithm in return.
          if (dynamic_cast<
                const parallel::distributed::Triangulation<dim, spacedim> *>(
                &dof_handler.get_triangulation()) != nullptr)
            for (typename hp::DoFHandler<dim, spacedim>::active_cell_iterator
                   cell = dof_handler.begin_active();
                 cell != dof_handler.end();
                 ++cell)
              if (cell->is_locally_owned() == false)
                for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
                     ++l)
                  cell->line(l)->set_user_flag();

          // An implementation of the algorithm described in the hp paper,
          // including the modification mentioned later in the "complications in
          // 3-d" subsections
          //
          // as explained there, we do something only if there are exactly 2
          // finite elements associated with an object. if there is only one,
          // then there is nothing to do anyway, and if there are 3 or more,
          // then we can get into trouble. note that this only happens for lines
          // in 3d and higher, and for quads only in 4d and higher, so this
          // isn't a particularly frequent case
          //
          // there is one case, however, that we would like to handle (see, for
          // example, the hp/crash_15 testcase): if we have
          // FESystem(FE_Q(2),FE_DGQ(i)) elements for a bunch of values 'i',
          // then we should be able to handle this because we can simply unify
          // *all* dofs, not only a some. so what we do is to first treat all
          // pairs of finite elements that have *identical* dofs, and then only
          // deal with those that are not identical of which we can handle at
          // most 2
          ::Table<2, std::unique_ptr<DoFIdentities>> line_dof_identities(
            dof_handler.fe_collection.size(), dof_handler.fe_collection.size());

          for (typename hp::DoFHandler<dim, spacedim>::active_cell_iterator
                 cell = dof_handler.begin_active();
               cell != dof_handler.end();
               ++cell)
            for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
              if (cell->line(l)->user_flag_set() == false)
                {
                  const typename hp::DoFHandler<dim, spacedim>::line_iterator
                    line = cell->line(l);
                  line->set_user_flag();

                  unsigned int unique_sets_of_dofs =
                    line->n_active_fe_indices();

                  // do a first loop over all sets of dofs and do identity
                  // uniquification
                  const unsigned int n_active_fe_indices =
                    line->n_active_fe_indices();
                  for (unsigned int f = 0; f < n_active_fe_indices; ++f)
                    for (unsigned int g = f + 1; g < n_active_fe_indices; ++g)
                      {
                        const unsigned int fe_index_1 =
                                             line->nth_active_fe_index(f),
                                           fe_index_2 =
                                             line->nth_active_fe_index(g);

                        if ((dof_handler.get_fe(fe_index_1).dofs_per_line ==
                             dof_handler.get_fe(fe_index_2).dofs_per_line) &&
                            (dof_handler.get_fe(fe_index_1).dofs_per_line > 0))
                          {
                            ensure_existence_of_dof_identities<1>(
                              dof_handler.get_fe(fe_index_1),
                              dof_handler.get_fe(fe_index_2),
                              line_dof_identities[fe_index_1][fe_index_2]);
                            // see if these sets of dofs are identical. the
                            // first condition for this is that indeed there are
                            // n identities
                            if (line_dof_identities[fe_index_1][fe_index_2]
                                  ->size() ==
                                dof_handler.get_fe(fe_index_1).dofs_per_line)
                              {
                                unsigned int i = 0;
                                for (; i < dof_handler.get_fe(fe_index_1)
                                             .dofs_per_line;
                                     ++i)
                                  if (((*(line_dof_identities[fe_index_1]
                                                             [fe_index_2]))[i]
                                         .first != i) &&
                                      ((*(line_dof_identities[fe_index_1]
                                                             [fe_index_2]))[i]
                                         .second != i))
                                    // not an identity
                                    break;

                                if (i == dof_handler.get_fe(fe_index_1)
                                           .dofs_per_line)
                                  {
                                    // The line dofs (i.e., the ones interior to
                                    // a line) of these two finite elements are
                                    // identical. Note that there could be
                                    // situations when one element still
                                    // dominates another, e.g.: FE_Q(2) x
                                    // FE_Nothing(dominate) vs FE_Q(2) x FE_Q(1)

                                    --unique_sets_of_dofs;

                                    for (unsigned int j = 0;
                                         j < dof_handler.get_fe(fe_index_1)
                                               .dofs_per_line;
                                         ++j)
                                      {
                                        const types::global_dof_index
                                          master_dof_index =
                                            line->dof_index(j, fe_index_1);
                                        const types::global_dof_index
                                          slave_dof_index =
                                            line->dof_index(j, fe_index_2);

                                        // if master dof was already
                                        // constrained, constrain to that one,
                                        // otherwise constrain slave to master
                                        if (dof_identities.find(
                                              master_dof_index) !=
                                            dof_identities.end())
                                          {
                                            Assert(dof_identities.find(
                                                     dof_identities
                                                       [master_dof_index]) ==
                                                     dof_identities.end(),
                                                   ExcInternalError());

                                            dof_identities[slave_dof_index] =
                                              dof_identities[master_dof_index];
                                          }
                                        else
                                          {
                                            Assert((dof_identities.find(
                                                      master_dof_index) ==
                                                    dof_identities.end()) ||
                                                     (dof_identities
                                                        [slave_dof_index] ==
                                                      master_dof_index),
                                                   ExcInternalError());

                                            dof_identities[slave_dof_index] =
                                              master_dof_index;
                                          }
                                      }
                                  }
                              }
                          }
                      }

                  // if at this point, there is only one unique set of dofs
                  // left, then we have taken care of everything above. if there
                  // are two, then we need to deal with them here. if there are
                  // more, then we punt, as described in the paper (and
                  // mentioned above)
                  // TODO: The check for 'dim==2' was inserted by intuition. It
                  // fixes
                  // the previous problems with @ref step_27 "step-27" in 3D. But an
                  // explanation for this is still required, and what we do here
                  // is not what we describe in the paper!.
                  if ((unique_sets_of_dofs == 2) && (dim == 2))
                    {
                      // find out which is the most dominating finite element of
                      // the ones that are used on this line
                      const unsigned int most_dominating_fe_index =
                        get_most_dominating_fe_index<dim, spacedim>(line);

                      // if we found the most dominating element, then use this
                      // to eliminate some of the degrees of freedom by
                      // identification. otherwise, the code that computes
                      // hanging node constraints will have to deal with it by
                      // computing appropriate constraints along this face/edge
                      if (most_dominating_fe_index !=
                          numbers::invalid_unsigned_int)
                        {
                          const unsigned int n_active_fe_indices =
                            line->n_active_fe_indices();

                          // loop over the indices of all the finite elements
                          // that are not dominating, and identify their dofs to
                          // the most dominating one
                          for (unsigned int f = 0; f < n_active_fe_indices; ++f)
                            if (line->nth_active_fe_index(f) !=
                                most_dominating_fe_index)
                              {
                                const unsigned int other_fe_index =
                                  line->nth_active_fe_index(f);

                                ensure_existence_of_dof_identities<1>(
                                  dof_handler.get_fe(most_dominating_fe_index),
                                  dof_handler.get_fe(other_fe_index),
                                  line_dof_identities[most_dominating_fe_index]
                                                     [other_fe_index]);

                                DoFIdentities &identities =
                                  *line_dof_identities[most_dominating_fe_index]
                                                      [other_fe_index];
                                for (unsigned int i = 0; i < identities.size();
                                     ++i)
                                  {
                                    const types::global_dof_index
                                      master_dof_index = line->dof_index(
                                        identities[i].first,
                                        most_dominating_fe_index);
                                    const types::global_dof_index
                                      slave_dof_index =
                                        line->dof_index(identities[i].second,
                                                        other_fe_index);

                                    Assert((dof_identities.find(
                                              master_dof_index) ==
                                            dof_identities.end()) ||
                                             (dof_identities[slave_dof_index] ==
                                              master_dof_index),
                                           ExcInternalError());

                                    dof_identities[slave_dof_index] =
                                      master_dof_index;
                                  }
                              }
                        }
                    }
                }

          // finally restore the user flags
          const_cast<::Triangulation<dim, spacedim> &>(
            dof_handler.get_triangulation())
            .load_user_flags_line(user_flags);

          return dof_identities;
        }



        template <int dim, int spacedim>
        static std::map<types::global_dof_index, types::global_dof_index>
        compute_quad_dof_identities(hp::DoFHandler<dim, spacedim> &)
        {
          // this function should only be called for dim<3 where there are
          // no quad dof identies. for dim>=3, the specialization below should
          // take care of it
          Assert(dim < 3, ExcInternalError());

          return std::map<types::global_dof_index, types::global_dof_index>();
        }


        template <int spacedim>
        static std::map<types::global_dof_index, types::global_dof_index>
          compute_quad_dof_identities(hp::DoFHandler<3, spacedim> &dof_handler)
        {
          const int dim = 3;

          std::map<types::global_dof_index, types::global_dof_index>
            dof_identities;


          // we will mark quads that we have already treated, so first
          // save and clear the user flags on quads and later restore
          // them
          std::vector<bool> user_flags;
          dof_handler.get_triangulation().save_user_flags_quad(user_flags);
          const_cast<::Triangulation<dim, spacedim> &>(
            dof_handler.get_triangulation())
            .clear_user_flags_quad();

          // exclude quads that bound cells we don't locally own, because
          // we do not have information about their dofs at this point.
          // this is, at least the case for
          // parallel::distributed::Triangulations.
          //
          // this means that we will not unify DoF indices between locally
          // owned cells and ghost cells, and this is different from what
          // we would do if the triangulation were not split into subdomains.
          // on the other hand, DoF unification is only an optimization: we
          // will still record these identities when we compute hanging
          // node constraints; we just end up with more DoFs than we would
          // if we unified DoF indices also between locally owned and ghost
          // cells, but we end up with a simpler algorithm in return.
          if (dynamic_cast<
                const parallel::distributed::Triangulation<dim, spacedim> *>(
                &dof_handler.get_triangulation()) != nullptr)
            for (typename hp::DoFHandler<dim, spacedim>::active_cell_iterator
                   cell = dof_handler.begin_active();
                 cell != dof_handler.end();
                 ++cell)
              if (cell->is_locally_owned() == false)
                for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell;
                     ++q)
                  cell->quad(q)->set_user_flag();


          // An implementation of the algorithm described in the hp
          // paper, including the modification mentioned later in the
          // "complications in 3-d" subsections
          //
          // as explained there, we do something only if there are
          // exactly 2 finite elements associated with an object. if
          // there is only one, then there is nothing to do anyway,
          // and if there are 3 or more, then we can get into
          // trouble. note that this only happens for lines in 3d and
          // higher, and for quads only in 4d and higher, so this
          // isn't a particularly frequent case
          ::Table<2, std::unique_ptr<DoFIdentities>> quad_dof_identities(
            dof_handler.fe_collection.size(), dof_handler.fe_collection.size());

          for (typename hp::DoFHandler<dim, spacedim>::active_cell_iterator
                 cell = dof_handler.begin_active();
               cell != dof_handler.end();
               ++cell)
            for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell; ++q)
              if ((cell->quad(q)->user_flag_set() == false) &&
                  (cell->quad(q)->n_active_fe_indices() == 2))
                {
                  const typename hp::DoFHandler<dim, spacedim>::quad_iterator
                    quad = cell->quad(q);
                  quad->set_user_flag();

                  // find out which is the most dominating finite
                  // element of the ones that are used on this quad
                  const unsigned int most_dominating_fe_index =
                    get_most_dominating_fe_index<dim, spacedim>(quad);

                  // if we found the most dominating element, then use
                  // this to eliminate some of the degrees of freedom
                  // by identification. otherwise, the code that
                  // computes hanging node constraints will have to
                  // deal with it by computing appropriate constraints
                  // along this face/edge
                  if (most_dominating_fe_index != numbers::invalid_unsigned_int)
                    {
                      const unsigned int n_active_fe_indices =
                        quad->n_active_fe_indices();

                      // loop over the indices of all the finite
                      // elements that are not dominating, and
                      // identify their dofs to the most dominating
                      // one
                      for (unsigned int f = 0; f < n_active_fe_indices; ++f)
                        if (quad->nth_active_fe_index(f) !=
                            most_dominating_fe_index)
                          {
                            const unsigned int other_fe_index =
                              quad->nth_active_fe_index(f);

                            ensure_existence_of_dof_identities<2>(
                              dof_handler.get_fe(most_dominating_fe_index),
                              dof_handler.get_fe(other_fe_index),
                              quad_dof_identities[most_dominating_fe_index]
                                                 [other_fe_index]);

                            DoFIdentities &identities =
                              *quad_dof_identities[most_dominating_fe_index]
                                                  [other_fe_index];
                            for (unsigned int i = 0; i < identities.size(); ++i)
                              {
                                const types::global_dof_index master_dof_index =
                                  quad->dof_index(identities[i].first,
                                                  most_dominating_fe_index);
                                const types::global_dof_index slave_dof_index =
                                  quad->dof_index(identities[i].second,
                                                  other_fe_index);

                                Assert((dof_identities.find(master_dof_index) ==
                                        dof_identities.end()) ||
                                         (dof_identities[slave_dof_index] ==
                                          master_dof_index),
                                       ExcInternalError());

                                dof_identities[slave_dof_index] =
                                  master_dof_index;
                              }
                          }
                    }
                }

          // finally restore the user flags
          const_cast<::Triangulation<dim, spacedim> &>(
            dof_handler.get_triangulation())
            .load_user_flags_quad(user_flags);

          return dof_identities;
        }



        template <int dim, int spacedim>
        static unsigned int
        unify_dof_indices(const DoFHandler<dim, spacedim> &,
                          const unsigned int n_dofs_before_identification,
                          const bool)
        {
          return n_dofs_before_identification;
        }



        template <int dim, int spacedim>
        static unsigned int
        unify_dof_indices(hp::DoFHandler<dim, spacedim> &dof_handler,
                          const unsigned int n_dofs_before_identification,
                          const bool         check_validity)
        {
          // compute the constraints that correspond to unifying
          // dof indices on vertices, lines, and quads. do so
          // in parallel
          std::map<types::global_dof_index, types::global_dof_index>
            all_constrained_indices[dim];

          {
            Threads::TaskGroup<> tasks;

            unsigned int i = 0;
            tasks += Threads::new_task([&, i]() {
              all_constrained_indices[i] =
                compute_vertex_dof_identities(dof_handler);
            });

            if (dim > 1)
              {
                ++i;
                tasks += Threads::new_task([&, i]() {
                  all_constrained_indices[i] =
                    compute_line_dof_identities(dof_handler);
                });
              }

            if (dim > 2)
              {
                ++i;
                tasks += Threads::new_task([&, i]() {
                  all_constrained_indices[i] =
                    compute_quad_dof_identities(dof_handler);
                });
              }

            tasks.join_all();
          }

          // create a vector that contains the new DoF indices; first preset the
          // ones that are identities as determined above, then enumerate the
          // rest
          std::vector<types::global_dof_index> new_dof_indices(
            n_dofs_before_identification, numbers::invalid_dof_index);

          for (const auto &constrained_dof_indices : all_constrained_indices)
            for (const auto &p : constrained_dof_indices)
              {
                Assert(new_dof_indices[p.first] == numbers::invalid_dof_index,
                       ExcInternalError());
                new_dof_indices[p.first] = p.second;
              }

          types::global_dof_index next_free_dof = 0;
          for (types::global_dof_index i = 0; i < n_dofs_before_identification;
               ++i)
            if (new_dof_indices[i] == numbers::invalid_dof_index)
              {
                new_dof_indices[i] = next_free_dof;
                ++next_free_dof;
              }

          // then loop over all those that are constrained and record the
          // new dof number for those:
          for (const auto &constrained_dof_indices : all_constrained_indices)
            for (const auto &p : constrained_dof_indices)
              {
                Assert(new_dof_indices[p.first] != numbers::invalid_dof_index,
                       ExcInternalError());

                new_dof_indices[p.first] = new_dof_indices[p.second];
              }

          for (types::global_dof_index i = 0; i < n_dofs_before_identification;
               ++i)
            {
              Assert(new_dof_indices[i] != numbers::invalid_dof_index,
                     ExcInternalError());
              Assert(new_dof_indices[i] < next_free_dof, ExcInternalError());
            }

          // finally, do the renumbering. verify that previous dof indices
          // were indeed all valid on all cells that we touch if we were
          // told to do so
          renumber_dofs(new_dof_indices,
                        IndexSet(0),
                        dof_handler,
                        check_validity);


          return next_free_dof;
        }



        template <class DoFHandlerType>
        static types::global_dof_index
        distribute_dofs(const types::subdomain_id subdomain_id,
                        DoFHandlerType &          dof_handler)
        {
          Assert(dof_handler.get_triangulation().n_levels() > 0,
                 ExcMessage("Empty triangulation"));

          // Step 1: distribute dofs on all cells, but definitely
          // exclude artificial cells
          types::global_dof_index next_free_dof = 0;
          typename DoFHandlerType::active_cell_iterator
            cell = dof_handler.begin_active(),
            endc = dof_handler.end();

          for (; cell != endc; ++cell)
            if (!cell->is_artificial())
              if ((subdomain_id == numbers::invalid_subdomain_id) ||
                  (cell->subdomain_id() == subdomain_id))
                next_free_dof =
                  Implementation::distribute_dofs_on_cell(dof_handler,
                                                          cell,
                                                          next_free_dof);

          // Step 2: unify dof indices in case this is an hp DoFHandler
          //
          // during unification, we need to renumber DoF indices. there,
          // we can check that all previous DoF indices were valid, but
          // this only makes sense if we really distributed DoFs on
          // all (non-artificial) cells above
          next_free_dof =
            unify_dof_indices(dof_handler,
                              next_free_dof,
                              /* check_validity = */
                              (subdomain_id == numbers::invalid_subdomain_id));

          update_all_active_cell_dof_indices_caches(dof_handler);

          return next_free_dof;
        }



        /* -------------- distribute_mg_dofs functionality ------------- */


        template <int dim, int spacedim>
        static types::global_dof_index
        distribute_mg_dofs_on_cell(
          const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell,
          types::global_dof_index next_free_dof,
          const std::integral_constant<int, 1> &)
        {
          // distribute dofs of vertices
          if (cell->get_fe().dofs_per_vertex > 0)
            for (unsigned int v = 0; v < GeometryInfo<1>::vertices_per_cell;
                 ++v)
              {
                typename DoFHandler<dim, spacedim>::level_cell_iterator
                  neighbor = cell->neighbor(v);

                if (neighbor.state() == IteratorState::valid)
                  {
                    // has neighbor already been processed?
                    if (neighbor->user_flag_set() &&
                        (neighbor->level() == cell->level()))
                      // copy dofs if the neighbor is on the same level (only
                      // then are mg dofs the same)
                      {
                        if (v == 0)
                          for (unsigned int d = 0;
                               d < cell->get_fe().dofs_per_vertex;
                               ++d)
                            cell->set_mg_vertex_dof_index(
                              cell->level(),
                              0,
                              d,
                              neighbor->mg_vertex_dof_index(cell->level(),
                                                            1,
                                                            d));
                        else
                          for (unsigned int d = 0;
                               d < cell->get_fe().dofs_per_vertex;
                               ++d)
                            cell->set_mg_vertex_dof_index(
                              cell->level(),
                              1,
                              d,
                              neighbor->mg_vertex_dof_index(cell->level(),
                                                            0,
                                                            d));

                        // next neighbor
                        continue;
                      }
                  }

                // otherwise: create dofs newly
                for (unsigned int d = 0; d < cell->get_fe().dofs_per_vertex;
                     ++d)
                  cell->set_mg_vertex_dof_index(cell->level(),
                                                v,
                                                d,
                                                next_free_dof++);
              }

          // dofs of line
          if (cell->get_fe().dofs_per_line > 0)
            for (unsigned int d = 0; d < cell->get_fe().dofs_per_line; ++d)
              cell->set_mg_dof_index(cell->level(), d, next_free_dof++);

          // note that this cell has been processed
          cell->set_user_flag();

          return next_free_dof;
        }



        template <int dim, int spacedim>
        static types::global_dof_index
        distribute_mg_dofs_on_cell(
          const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell,
          types::global_dof_index next_free_dof,
          const std::integral_constant<int, 2> &)
        {
          if (cell->get_fe().dofs_per_vertex > 0)
            // number dofs on vertices
            for (unsigned int vertex = 0;
                 vertex < GeometryInfo<2>::vertices_per_cell;
                 ++vertex)
              // check whether dofs for this
              // vertex have been distributed
              // (only check the first dof)
              if (cell->mg_vertex_dof_index(cell->level(), vertex, 0) ==
                  numbers::invalid_dof_index)
                for (unsigned int d = 0; d < cell->get_fe().dofs_per_vertex;
                     ++d)
                  cell->set_mg_vertex_dof_index(cell->level(),
                                                vertex,
                                                d,
                                                next_free_dof++);

          // for the four sides
          if (cell->get_fe().dofs_per_line > 0)
            for (unsigned int side = 0; side < GeometryInfo<2>::faces_per_cell;
                 ++side)
              {
                typename DoFHandler<dim, spacedim>::line_iterator line =
                  cell->line(side);

                // distribute dofs if necessary: check whether line dof is
                // already numbered (check only first dof)
                if (line->mg_dof_index(cell->level(), 0) ==
                    numbers::invalid_dof_index)
                  // if not: distribute dofs
                  for (unsigned int d = 0; d < cell->get_fe().dofs_per_line;
                       ++d)
                    line->set_mg_dof_index(cell->level(), d, next_free_dof++);
              }


          // dofs of quad
          if (cell->get_fe().dofs_per_quad > 0)
            for (unsigned int d = 0; d < cell->get_fe().dofs_per_quad; ++d)
              cell->set_mg_dof_index(cell->level(), d, next_free_dof++);


          // note that this cell has been processed
          cell->set_user_flag();

          return next_free_dof;
        }



        template <int dim, int spacedim>
        static types::global_dof_index
        distribute_mg_dofs_on_cell(
          const typename DoFHandler<dim, spacedim>::level_cell_iterator &cell,
          types::global_dof_index next_free_dof,
          const std::integral_constant<int, 3> &)
        {
          if (cell->get_fe().dofs_per_vertex > 0)
            // number dofs on vertices
            for (unsigned int vertex = 0;
                 vertex < GeometryInfo<3>::vertices_per_cell;
                 ++vertex)
              // check whether dofs for this vertex have been distributed
              // (only check the first dof)
              if (cell->mg_vertex_dof_index(cell->level(), vertex, 0) ==
                  numbers::invalid_dof_index)
                for (unsigned int d = 0; d < cell->get_fe().dofs_per_vertex;
                     ++d)
                  cell->set_mg_vertex_dof_index(cell->level(),
                                                vertex,
                                                d,
                                                next_free_dof++);

          // for the lines
          if (cell->get_fe().dofs_per_line > 0)
            for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell; ++l)
              {
                typename DoFHandler<dim, spacedim>::line_iterator line =
                  cell->line(l);

                // distribute dofs if necessary:
                // check whether line dof is already
                // numbered (check only first dof)
                if (line->mg_dof_index(cell->level(), 0) ==
                    numbers::invalid_dof_index)
                  // if not: distribute dofs
                  for (unsigned int d = 0; d < cell->get_fe().dofs_per_line;
                       ++d)
                    line->set_mg_dof_index(cell->level(), d, next_free_dof++);
              }

          // for the quads
          if (cell->get_fe().dofs_per_quad > 0)
            for (unsigned int q = 0; q < GeometryInfo<3>::quads_per_cell; ++q)
              {
                typename DoFHandler<dim, spacedim>::quad_iterator quad =
                  cell->quad(q);

                // distribute dofs if necessary:
                // check whether line dof is already
                // numbered (check only first dof)
                if (quad->mg_dof_index(cell->level(), 0) ==
                    numbers::invalid_dof_index)
                  // if not: distribute dofs
                  for (unsigned int d = 0; d < cell->get_fe().dofs_per_quad;
                       ++d)
                    quad->set_mg_dof_index(cell->level(), d, next_free_dof++);
              }


          // dofs of cell
          if (cell->get_fe().dofs_per_hex > 0)
            for (unsigned int d = 0; d < cell->get_fe().dofs_per_hex; ++d)
              cell->set_mg_dof_index(cell->level(), d, next_free_dof++);


          // note that this cell has been processed
          cell->set_user_flag();

          return next_free_dof;
        }



        // same for the hp::DoFHandler
        template <int spacedim>
        static types::global_dof_index
        distribute_mg_dofs_on_cell(
          const hp::DoFHandler<1, spacedim> &dof_handler,
          const typename hp::DoFHandler<1, spacedim>::active_cell_iterator
            &                     cell,
          types::global_dof_index next_free_dof)
        {
          (void)dof_handler;
          (void)cell;
          (void)next_free_dof;
          return 0;
        }



        template <int spacedim>
        static types::global_dof_index
        distribute_mg_dofs_on_cell(
          const hp::DoFHandler<2, spacedim> &dof_handler,
          const typename hp::DoFHandler<2, spacedim>::active_cell_iterator
            &                     cell,
          types::global_dof_index next_free_dof)
        {
          (void)dof_handler;
          (void)cell;
          (void)next_free_dof;
          return 0;
        }



        template <int spacedim>
        static types::global_dof_index
        distribute_mg_dofs_on_cell(
          const hp::DoFHandler<3, spacedim> &dof_handler,
          const typename hp::DoFHandler<3, spacedim>::active_cell_iterator
            &                     cell,
          types::global_dof_index next_free_dof)
        {
          (void)dof_handler;
          (void)cell;
          (void)next_free_dof;
          return 0;
        }



        template <class DoFHandlerType>
        static types::global_dof_index
        distribute_dofs_on_level(const types::subdomain_id level_subdomain_id,
                                 DoFHandlerType &          dof_handler,
                                 const unsigned int        level)
        {
          const unsigned int dim      = DoFHandlerType::dimension;
          const unsigned int spacedim = DoFHandlerType::space_dimension;

          const ::Triangulation<dim, spacedim> &tria =
            dof_handler.get_triangulation();
          Assert(tria.n_levels() > 0, ExcMessage("Empty triangulation"));
          if (level >= tria.n_levels())
            return 0; // this is allowed for multigrid

          // Clear user flags because we will need them. But first we save
          // them and make sure that we restore them later such that at
          // the end of this function the Triangulation will be in the
          // same state as it was at the beginning of this function.
          std::vector<bool> user_flags;
          tria.save_user_flags(user_flags);
          const_cast<::Triangulation<dim, spacedim> &>(tria)
            .clear_user_flags();

          types::global_dof_index next_free_dof = 0;
          typename DoFHandler<dim, spacedim>::level_cell_iterator
            cell = dof_handler.begin(level),
            endc = dof_handler.end(level);

          for (; cell != endc; ++cell)
            if ((level_subdomain_id == numbers::invalid_subdomain_id) ||
                (cell->level_subdomain_id() == level_subdomain_id))
              next_free_dof =
                Implementation::distribute_mg_dofs_on_cell<dim, spacedim>(
                  cell, next_free_dof, std::integral_constant<int, dim>());

          // finally restore the user flags
          const_cast<::Triangulation<dim, spacedim> &>(tria)
            .load_user_flags(user_flags);

          return next_free_dof;
        }



        /* --------------------- renumber_dofs functionality ---------------- */


        template <int dim, int spacedim>
        static void
        renumber_vertex_dofs(
          const std::vector<types::global_dof_index> &new_numbers,
          const IndexSet &                            indices,
          DoFHandler<dim, spacedim> &                 dof_handler,
          const bool                                  check_validity)
        {
          // we can not use cell iterators in this function since then
          // we would renumber the dofs on the interface of two cells
          // more than once. Anyway, this way it's not only more
          // correct but also faster; note, however, that dof numbers
          // may be invalid_dof_index, namely when the appropriate
          // vertex/line/etc is unused
          for (std::vector<types::global_dof_index>::iterator i =
                 dof_handler.vertex_dofs.begin();
               i != dof_handler.vertex_dofs.end();
               ++i)
            if (*i != numbers::invalid_dof_index)
              *i = (indices.size() == 0) ?
                     (new_numbers[*i]) :
                     (new_numbers[indices.index_within_set(*i)]);
            else if (check_validity)
              // if index is invalid_dof_index: check if this one
              // really is unused
              Assert(dof_handler.get_triangulation().vertex_used(
                       (i - dof_handler.vertex_dofs.begin()) /
                       dof_handler.get_fe().dofs_per_vertex) == false,
                     ExcInternalError());
        }



        template <int dim, int spacedim>
        static void
        renumber_cell_dofs(
          const std::vector<types::global_dof_index> &new_numbers,
          const IndexSet &                            indices,
          DoFHandler<dim, spacedim> &                 dof_handler)
        {
          for (unsigned int level = 0; level < dof_handler.levels.size();
               ++level)
            for (std::vector<types::global_dof_index>::iterator i =
                   dof_handler.levels[level]->dof_object.dofs.begin();
                 i != dof_handler.levels[level]->dof_object.dofs.end();
                 ++i)
              if (*i != numbers::invalid_dof_index)
                *i = ((indices.size() == 0) ?
                        new_numbers[*i] :
                        new_numbers[indices.index_within_set(*i)]);
        }



        template <int spacedim>
        static void
        renumber_face_dofs(
          const std::vector<types::global_dof_index> & /*new_numbers*/,
          const IndexSet & /*indices*/,
          DoFHandler<1, spacedim> & /*dof_handler*/)
        {
          // nothing to do in 1d since there are no separate faces
        }



        template <int spacedim>
        static void
        renumber_face_dofs(
          const std::vector<types::global_dof_index> &new_numbers,
          const IndexSet &                            indices,
          DoFHandler<2, spacedim> &                   dof_handler)
        {
          // treat dofs on lines
          for (std::vector<types::global_dof_index>::iterator i =
                 dof_handler.faces->lines.dofs.begin();
               i != dof_handler.faces->lines.dofs.end();
               ++i)
            if (*i != numbers::invalid_dof_index)
              *i = ((indices.size() == 0) ?
                      new_numbers[*i] :
                      new_numbers[indices.index_within_set(*i)]);
        }



        template <int spacedim>
        static void
        renumber_face_dofs(
          const std::vector<types::global_dof_index> &new_numbers,
          const IndexSet &                            indices,
          DoFHandler<3, spacedim> &                   dof_handler)
        {
          // treat dofs on lines
          for (std::vector<types::global_dof_index>::iterator i =
                 dof_handler.faces->lines.dofs.begin();
               i != dof_handler.faces->lines.dofs.end();
               ++i)
            if (*i != numbers::invalid_dof_index)
              *i = ((indices.size() == 0) ?
                      new_numbers[*i] :
                      new_numbers[indices.index_within_set(*i)]);

          // treat dofs on quads
          for (std::vector<types::global_dof_index>::iterator i =
                 dof_handler.faces->quads.dofs.begin();
               i != dof_handler.faces->quads.dofs.end();
               ++i)
            if (*i != numbers::invalid_dof_index)
              *i = ((indices.size() == 0) ?
                      new_numbers[*i] :
                      new_numbers[indices.index_within_set(*i)]);
        }



        template <int dim, int spacedim>
        static void
        renumber_vertex_dofs(
          const std::vector<types::global_dof_index> &new_numbers,
          const IndexSet &                            indices,
          hp::DoFHandler<dim, spacedim> &             dof_handler,
          const bool                                  check_validity)
        {
          for (unsigned int vertex_index = 0;
               vertex_index < dof_handler.get_triangulation().n_vertices();
               ++vertex_index)
            {
              const unsigned int n_active_fe_indices =
                ::internal::DoFAccessorImplementation::Implementation::
                  n_active_vertex_fe_indices(dof_handler, vertex_index);

              // if this vertex is unused, then we really ought not to have
              // allocated any space for it, i.e., n_active_fe_indices should be
              // zero, and there is no space to actually store dof indices for
              // this vertex
              if (dof_handler.get_triangulation().vertex_used(vertex_index) ==
                  false)
                Assert(n_active_fe_indices == 0, ExcInternalError());

              // otherwise the vertex is used; it may still not hold any dof
              // indices if it is located on an artificial cell and not adjacent
              // to a ghost cell, but in that case there is simply nothing for
              // us to do
              for (unsigned int f = 0; f < n_active_fe_indices; ++f)
                {
                  const unsigned int fe_index =
                    ::internal::DoFAccessorImplementation::
                      Implementation::nth_active_vertex_fe_index(dof_handler,
                                                                 vertex_index,
                                                                 f);

                  for (unsigned int d = 0;
                       d < dof_handler.get_fe(fe_index).dofs_per_vertex;
                       ++d)
                    {
                      const types::global_dof_index old_dof_index =
                        ::internal::DoFAccessorImplementation::
                          Implementation::get_vertex_dof_index(dof_handler,
                                                               vertex_index,
                                                               fe_index,
                                                               d);

                      // if check_validity was set, then we are to verify that
                      // the previous indices were all valid. this really should
                      // be the case: we allocated space for these vertex dofs,
                      // i.e., at least one adjacent cell has a valid
                      // active_fe_index, so there are DoFs that really live
                      // on this vertex. if check_validity is set, then we
                      // must make sure that they have been set to something
                      // useful
                      if (check_validity)
                        Assert(old_dof_index != numbers::invalid_dof_index,
                               ExcInternalError());

                      if (old_dof_index != numbers::invalid_dof_index)
                        ::internal::DoFAccessorImplementation::
                          Implementation::set_vertex_dof_index(
                            dof_handler,
                            vertex_index,
                            fe_index,
                            d,
                            (indices.size() == 0) ?
                              (new_numbers[old_dof_index]) :
                              (new_numbers[indices.index_within_set(
                                old_dof_index)]));
                    }
                }
            }
        }



        template <int dim, int spacedim>
        static void
        renumber_cell_dofs(
          const std::vector<types::global_dof_index> &new_numbers,
          const IndexSet &                            indices,
          hp::DoFHandler<dim, spacedim> &             dof_handler)
        {
          for (typename hp::DoFHandler<dim, spacedim>::active_cell_iterator
                 cell = dof_handler.begin_active();
               cell != dof_handler.end();
               ++cell)
            if (!cell->is_artificial())
              {
                const unsigned int fe_index = cell->active_fe_index();

                for (unsigned int d = 0;
                     d < dof_handler.get_fe(fe_index)
                           .template n_dofs_per_object<dim>();
                     ++d)
                  {
                    const types::global_dof_index old_dof_index =
                      cell->dof_index(d, fe_index);
                    if (old_dof_index != numbers::invalid_dof_index)
                      cell->set_dof_index(
                        d,
                        (indices.size() == 0) ?
                          (new_numbers[old_dof_index]) :
                          (new_numbers[indices.index_within_set(
                            old_dof_index)]),
                        fe_index);
                  }
              }
        }



        template <int spacedim>
        static void
        renumber_face_dofs(
          const std::vector<types::global_dof_index> & /*new_numbers*/,
          const IndexSet & /*indices*/,
          hp::DoFHandler<1, spacedim> & /*dof_handler*/)
        {
          // nothing to do in 1d since there are no separate faces -- we've
          // already taken care of this when dealing with the vertices
        }



        template <int spacedim>
        static void
        renumber_face_dofs(
          const std::vector<types::global_dof_index> &new_numbers,
          const IndexSet &                            indices,
          hp::DoFHandler<2, spacedim> &               dof_handler)
        {
          const unsigned int dim = 2;

          // deal with DoFs on lines
          {
            // save user flags on lines so we can use them to mark lines
            // we've already treated
            std::vector<bool> saved_line_user_flags;
            const_cast<::Triangulation<dim, spacedim> &>(
              dof_handler.get_triangulation())
              .save_user_flags_line(saved_line_user_flags);
            const_cast<::Triangulation<dim, spacedim> &>(
              dof_handler.get_triangulation())
              .clear_user_flags_line();

            for (typename hp::DoFHandler<dim, spacedim>::active_cell_iterator
                   cell = dof_handler.begin_active();
                 cell != dof_handler.end();
                 ++cell)
              if (!cell->is_artificial())
                for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
                     ++l)
                  if (cell->line(l)->user_flag_set() == false)
                    {
                      const typename hp::DoFHandler<dim,
                                                    spacedim>::line_iterator
                        line = cell->line(l);
                      line->set_user_flag();

                      const unsigned int n_active_fe_indices =
                        line->n_active_fe_indices();

                      for (unsigned int f = 0; f < n_active_fe_indices; ++f)
                        {
                          const unsigned int fe_index =
                            line->nth_active_fe_index(f);

                          for (unsigned int d = 0;
                               d < dof_handler.get_fe(fe_index).dofs_per_line;
                               ++d)
                            {
                              const types::global_dof_index old_dof_index =
                                line->dof_index(d, fe_index);
                              if (old_dof_index != numbers::invalid_dof_index)
                                line->set_dof_index(
                                  d,
                                  (indices.size() == 0) ?
                                    (new_numbers[old_dof_index]) :
                                    (new_numbers[indices.index_within_set(
                                      old_dof_index)]),
                                  fe_index);
                            }
                        }
                    }

            // at the end, restore the user
            // flags for the lines
            const_cast<::Triangulation<dim, spacedim> &>(
              dof_handler.get_triangulation())
              .load_user_flags_line(saved_line_user_flags);
          }
        }



        template <int spacedim>
        static void
        renumber_face_dofs(
          const std::vector<types::global_dof_index> &new_numbers,
          const IndexSet &                            indices,
          hp::DoFHandler<3, spacedim> &               dof_handler)
        {
          const unsigned int dim = 3;

          // deal with DoFs on lines
          {
            // save user flags on lines so we can use them to mark lines
            // we've already treated
            std::vector<bool> saved_line_user_flags;
            const_cast<::Triangulation<dim, spacedim> &>(
              dof_handler.get_triangulation())
              .save_user_flags_line(saved_line_user_flags);
            const_cast<::Triangulation<dim, spacedim> &>(
              dof_handler.get_triangulation())
              .clear_user_flags_line();

            for (typename hp::DoFHandler<dim, spacedim>::active_cell_iterator
                   cell = dof_handler.begin_active();
                 cell != dof_handler.end();
                 ++cell)
              if (!cell->is_artificial())
                for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell;
                     ++l)
                  if (cell->line(l)->user_flag_set() == false)
                    {
                      const typename hp::DoFHandler<dim,
                                                    spacedim>::line_iterator
                        line = cell->line(l);
                      line->set_user_flag();

                      const unsigned int n_active_fe_indices =
                        line->n_active_fe_indices();

                      for (unsigned int f = 0; f < n_active_fe_indices; ++f)
                        {
                          const unsigned int fe_index =
                            line->nth_active_fe_index(f);

                          for (unsigned int d = 0;
                               d < dof_handler.get_fe(fe_index).dofs_per_line;
                               ++d)
                            {
                              const types::global_dof_index old_dof_index =
                                line->dof_index(d, fe_index);
                              if (old_dof_index != numbers::invalid_dof_index)
                                line->set_dof_index(
                                  d,
                                  (indices.size() == 0) ?
                                    (new_numbers[old_dof_index]) :
                                    (new_numbers[indices.index_within_set(
                                      old_dof_index)]),
                                  fe_index);
                            }
                        }
                    }

            // at the end, restore the user
            // flags for the lines
            const_cast<::Triangulation<dim, spacedim> &>(
              dof_handler.get_triangulation())
              .load_user_flags_line(saved_line_user_flags);
          }

          // then deal with dofs on quads
          {
            std::vector<bool> saved_quad_user_flags;
            const_cast<::Triangulation<dim, spacedim> &>(
              dof_handler.get_triangulation())
              .save_user_flags_quad(saved_quad_user_flags);
            const_cast<::Triangulation<dim, spacedim> &>(
              dof_handler.get_triangulation())
              .clear_user_flags_quad();

            for (typename hp::DoFHandler<dim, spacedim>::active_cell_iterator
                   cell = dof_handler.begin_active();
                 cell != dof_handler.end();
                 ++cell)
              if (!cell->is_artificial())
                for (unsigned int q = 0; q < GeometryInfo<dim>::quads_per_cell;
                     ++q)
                  if (cell->quad(q)->user_flag_set() == false)
                    {
                      const typename hp::DoFHandler<dim,
                                                    spacedim>::quad_iterator
                        quad = cell->quad(q);
                      quad->set_user_flag();

                      const unsigned int n_active_fe_indices =
                        quad->n_active_fe_indices();

                      for (unsigned int f = 0; f < n_active_fe_indices; ++f)
                        {
                          const unsigned int fe_index =
                            quad->nth_active_fe_index(f);

                          for (unsigned int d = 0;
                               d < dof_handler.get_fe(fe_index).dofs_per_quad;
                               ++d)
                            {
                              const types::global_dof_index old_dof_index =
                                quad->dof_index(d, fe_index);
                              if (old_dof_index != numbers::invalid_dof_index)
                                quad->set_dof_index(
                                  d,
                                  (indices.size() == 0) ?
                                    (new_numbers[old_dof_index]) :
                                    (new_numbers[indices.index_within_set(
                                      old_dof_index)]),
                                  fe_index);
                            }
                        }
                    }

            // at the end, restore the user flags for the quads
            const_cast<::Triangulation<dim, spacedim> &>(
              dof_handler.get_triangulation())
              .load_user_flags_quad(saved_quad_user_flags);
          }
        }



        template <class DoFHandlerType>
        static void
        renumber_dofs(const std::vector<types::global_dof_index> &new_numbers,
                      const IndexSet &                            indices,
                      DoFHandlerType &                            dof_handler,
                      const bool check_validity)
        {
          if (DoFHandlerType::dimension == 1)
            Assert(indices == IndexSet(0), ExcNotImplemented());

          // renumber DoF indices on vertices, cells, and faces. this
          // can be done in parallel because the respective functions
          // work on separate data structures
          Threads::TaskGroup<> tasks;
          tasks += Threads::new_task([&]() {
            renumber_vertex_dofs(new_numbers,
                                 indices,
                                 dof_handler,
                                 check_validity);
          });
          tasks += Threads::new_task(
            [&]() { renumber_face_dofs(new_numbers, indices, dof_handler); });
          tasks += Threads::new_task(
            [&]() { renumber_cell_dofs(new_numbers, indices, dof_handler); });
          tasks.join_all();

          // update the cache used for cell dof indices
          update_all_active_cell_dof_indices_caches(dof_handler);
        }



        /* --------------------- renumber_mg_dofs functionality ----------------
         */

        template <int dim, int spacedim>
        static void
        renumber_vertex_mg_dofs(
          const std::vector<::types::global_dof_index> &new_numbers,
          const IndexSet &                                    indices,
          DoFHandler<dim, spacedim> &                         dof_handler,
          const unsigned int                                  level,
          const bool                                          check_validity)
        {
          Assert(level < dof_handler.get_triangulation().n_levels(),
                 ExcInternalError());

          for (typename std::vector<
                 typename DoFHandler<dim, spacedim>::MGVertexDoFs>::iterator i =
                 dof_handler.mg_vertex_dofs.begin();
               i != dof_handler.mg_vertex_dofs.end();
               ++i)
            // if the present vertex lives on the current level
            if ((i->get_coarsest_level() <= level) &&
                (i->get_finest_level() >= level))
              for (unsigned int d = 0; d < dof_handler.get_fe().dofs_per_vertex;
                   ++d)
                {
                  const ::types::global_dof_index idx =
                    i->get_index(level,
                                 d,
                                 dof_handler.get_fe().dofs_per_vertex);

                  if (check_validity)
                    Assert(idx != numbers::invalid_dof_index,
                           ExcInternalError());

                  if (idx != numbers::invalid_dof_index)
                    i->set_index(
                      level,
                      d,
                      dof_handler.get_fe().dofs_per_vertex,
                      (indices.size() == 0) ?
                        (new_numbers[idx]) :
                        (new_numbers[indices.index_within_set(idx)]));
                }
        }



        template <int dim, int spacedim>
        static void
        renumber_cell_mg_dofs(
          const std::vector<::types::global_dof_index> &new_numbers,
          const IndexSet &                                    indices,
          DoFHandler<dim, spacedim> &                         dof_handler,
          const unsigned int                                  level)
        {
          for (std::vector<types::global_dof_index>::iterator i =
                 dof_handler.mg_levels[level]->dof_object.dofs.begin();
               i != dof_handler.mg_levels[level]->dof_object.dofs.end();
               ++i)
            {
              if (*i != numbers::invalid_dof_index)
                {
                  Assert((indices.size() > 0 ? indices.is_element(*i) :
                                               (*i < new_numbers.size())),
                         ExcInternalError());
                  *i = (indices.size() == 0) ?
                         (new_numbers[*i]) :
                         (new_numbers[indices.index_within_set(*i)]);
                }
            }
        }



        template <int spacedim>
        static void
        renumber_face_mg_dofs(
          const std::vector<types::global_dof_index> & /*new_numbers*/,
          const IndexSet & /*indices*/,
          DoFHandler<1, spacedim> & /*dof_handler*/,
          const unsigned int /*level*/,
          const bool /*check_validity*/)
        {
          // nothing to do in 1d because there are no separate faces
        }



        template <int spacedim>
        static void
        renumber_face_mg_dofs(
          const std::vector<::types::global_dof_index> &new_numbers,
          const IndexSet &                                    indices,
          DoFHandler<2, spacedim> &                           dof_handler,
          const unsigned int                                  level,
          const bool                                          check_validity)
        {
          if (dof_handler.get_fe().dofs_per_line > 0)
            {
              // save user flags as they will be modified
              std::vector<bool> user_flags;
              dof_handler.get_triangulation().save_user_flags(user_flags);
              const_cast<::Triangulation<2, spacedim> &>(
                dof_handler.get_triangulation())
                .clear_user_flags();

              // flag all lines adjacent to cells of the current
              // level, as those lines logically belong to the same
              // level as the cell, at least for for isotropic
              // refinement
              typename DoFHandler<2, spacedim>::level_cell_iterator cell,
                endc = dof_handler.end(level);
              for (cell = dof_handler.begin(level); cell != endc; ++cell)
                for (unsigned int line = 0;
                     line < GeometryInfo<2>::faces_per_cell;
                     ++line)
                  cell->face(line)->set_user_flag();

              for (typename DoFHandler<2, spacedim>::cell_iterator cell =
                     dof_handler.begin();
                   cell != dof_handler.end();
                   ++cell)
                for (unsigned int l = 0; l < GeometryInfo<2>::lines_per_cell;
                     ++l)
                  if (cell->line(l)->user_flag_set())
                    {
                      for (unsigned int d = 0;
                           d < dof_handler.get_fe().dofs_per_line;
                           ++d)
                        {
                          const ::types::global_dof_index idx =
                            cell->line(l)->mg_dof_index(level, d);
                          if (check_validity)
                            Assert(idx != numbers::invalid_dof_index,
                                   ExcInternalError());

                          if (idx != numbers::invalid_dof_index)
                            cell->line(l)->set_mg_dof_index(
                              level,
                              d,
                              ((indices.size() == 0) ?
                                 new_numbers[idx] :
                                 new_numbers[indices.index_within_set(idx)]));
                        }
                      cell->line(l)->clear_user_flag();
                    }
              // finally, restore user flags
              const_cast<::Triangulation<2, spacedim> &>(
                dof_handler.get_triangulation())
                .load_user_flags(user_flags);
            }
        }



        template <int spacedim>
        static void
        renumber_face_mg_dofs(
          const std::vector<::types::global_dof_index> &new_numbers,
          const IndexSet &                                    indices,
          DoFHandler<3, spacedim> &                           dof_handler,
          const unsigned int                                  level,
          const bool                                          check_validity)
        {
          if (dof_handler.get_fe().dofs_per_line > 0 ||
              dof_handler.get_fe().dofs_per_quad > 0)
            {
              // save user flags as they will be modified
              std::vector<bool> user_flags;
              dof_handler.get_triangulation().save_user_flags(user_flags);
              const_cast<::Triangulation<3, spacedim> &>(
                dof_handler.get_triangulation())
                .clear_user_flags();

              // flag all lines adjacent to cells of the current
              // level, as those lines logically belong to the same
              // level as the cell, at least for isotropic refinement
              typename DoFHandler<3, spacedim>::level_cell_iterator cell,
                endc = dof_handler.end(level);
              for (cell = dof_handler.begin(level); cell != endc; ++cell)
                for (unsigned int line = 0;
                     line < GeometryInfo<3>::lines_per_cell;
                     ++line)
                  cell->line(line)->set_user_flag();

              for (typename DoFHandler<3, spacedim>::cell_iterator cell =
                     dof_handler.begin();
                   cell != dof_handler.end();
                   ++cell)
                for (unsigned int l = 0; l < GeometryInfo<3>::lines_per_cell;
                     ++l)
                  if (cell->line(l)->user_flag_set())
                    {
                      for (unsigned int d = 0;
                           d < dof_handler.get_fe().dofs_per_line;
                           ++d)
                        {
                          const ::types::global_dof_index idx =
                            cell->line(l)->mg_dof_index(level, d);
                          if (check_validity)
                            Assert(idx != numbers::invalid_dof_index,
                                   ExcInternalError());

                          if (idx != numbers::invalid_dof_index)
                            cell->line(l)->set_mg_dof_index(
                              level,
                              d,
                              ((indices.size() == 0) ?
                                 new_numbers[idx] :
                                 new_numbers[indices.index_within_set(idx)]));
                        }
                      cell->line(l)->clear_user_flag();
                    }

              // flag all quads adjacent to cells of the current level, as
              // those quads logically belong to the same level as the cell,
              // at least for isotropic refinement
              for (cell = dof_handler.begin(level); cell != endc; ++cell)
                for (unsigned int quad = 0;
                     quad < GeometryInfo<3>::quads_per_cell;
                     ++quad)
                  cell->quad(quad)->set_user_flag();

              for (typename DoFHandler<3, spacedim>::cell_iterator cell =
                     dof_handler.begin();
                   cell != dof_handler.end();
                   ++cell)
                for (unsigned int l = 0; l < GeometryInfo<3>::quads_per_cell;
                     ++l)
                  if (cell->quad(l)->user_flag_set())
                    {
                      for (unsigned int d = 0;
                           d < dof_handler.get_fe().dofs_per_quad;
                           ++d)
                        {
                          const ::types::global_dof_index idx =
                            cell->quad(l)->mg_dof_index(level, d);
                          if (check_validity)
                            Assert(idx != numbers::invalid_dof_index,
                                   ExcInternalError());

                          if (idx != numbers::invalid_dof_index)
                            cell->quad(l)->set_mg_dof_index(
                              level,
                              d,
                              ((indices.size() == 0) ?
                                 new_numbers[idx] :
                                 new_numbers[indices.index_within_set(idx)]));
                        }
                      cell->quad(l)->clear_user_flag();
                    }

              // finally, restore user flags
              const_cast<::Triangulation<3, spacedim> &>(
                dof_handler.get_triangulation())
                .load_user_flags(user_flags);
            }
        }



        template <int dim, int spacedim>
        static void
        renumber_mg_dofs(
          const std::vector<::types::global_dof_index> &new_numbers,
          const IndexSet &                                    indices,
          DoFHandler<dim, spacedim> &                         dof_handler,
          const unsigned int                                  level,
          const bool                                          check_validity)
        {
          Assert(level < dof_handler.get_triangulation().n_global_levels(),
                 ExcInternalError());

          // renumber DoF indices on vertices, cells, and faces. this
          // can be done in parallel because the respective functions
          // work on separate data structures
          Threads::TaskGroup<> tasks;
          tasks += Threads::new_task([&]() {
            renumber_vertex_mg_dofs(
              new_numbers, indices, dof_handler, level, check_validity);
          });
          tasks += Threads::new_task([&]() {
            renumber_face_mg_dofs(
              new_numbers, indices, dof_handler, level, check_validity);
          });
          tasks += Threads::new_task([&]() {
            renumber_cell_mg_dofs(new_numbers, indices, dof_handler, level);
          });
          tasks.join_all();
        }



        template <int dim, int spacedim>
        static void
        renumber_mg_dofs(
          const std::vector<::types::global_dof_index> & /*new_numbers*/,
          const IndexSet & /*indices*/,
          hp::DoFHandler<dim, spacedim> & /*dof_handler*/,
          const unsigned int /*level*/,
          const bool /*check_validity*/)
        {
          Assert(false, ExcNotImplemented());
        }
      };



      /* --------------------- class Sequential ---------------- */



      template <class DoFHandlerType>
      Sequential<DoFHandlerType>::Sequential(DoFHandlerType &dof_handler)
        : dof_handler(&dof_handler)
      {}



      template <class DoFHandlerType>
      NumberCache
      Sequential<DoFHandlerType>::distribute_dofs() const
      {
        const types::global_dof_index n_dofs =
          Implementation::distribute_dofs(numbers::invalid_subdomain_id,
                                          *dof_handler);

        // return a sequential, complete index set
        return NumberCache(n_dofs);
      }



      template <class DoFHandlerType>
      std::vector<NumberCache>
      Sequential<DoFHandlerType>::distribute_mg_dofs() const
      {
        std::vector<bool> user_flags;
        dof_handler->get_triangulation().save_user_flags(user_flags);

        const_cast<::Triangulation<DoFHandlerType::dimension,
                                         DoFHandlerType::space_dimension> &>(
          dof_handler->get_triangulation())
          .clear_user_flags();

        std::vector<NumberCache> number_caches;
        number_caches.reserve(dof_handler->get_triangulation().n_levels());
        for (unsigned int level = 0;
             level < dof_handler->get_triangulation().n_levels();
             ++level)
          {
            // first distribute dofs on this level
            const types::global_dof_index n_level_dofs =
              Implementation::distribute_dofs_on_level(
                numbers::invalid_subdomain_id, *dof_handler, level);

            // then add a complete, sequential index set
            number_caches.emplace_back(NumberCache(n_level_dofs));
          }

        const_cast<::Triangulation<DoFHandlerType::dimension,
                                         DoFHandlerType::space_dimension> &>(
          dof_handler->get_triangulation())
          .load_user_flags(user_flags);

        return number_caches;
      }



      template <class DoFHandlerType>
      NumberCache
      Sequential<DoFHandlerType>::renumber_dofs(
        const std::vector<types::global_dof_index> &new_numbers) const
      {
        Implementation::renumber_dofs(new_numbers,
                                      IndexSet(0),
                                      *dof_handler,
                                      true);

        // return a sequential, complete index set. take into account that the
        // number of DoF indices may in fact be smaller than there were before
        // if some previously separately numbered dofs have been identified.
        // this is, for example, what the hp::DoFHandler does: it first
        // enumerates all DoFs on cells independently, and then unifies
        // some located at vertices or faces; this leaves us with fewer
        // DoFs than there were before, so use the largest index as
        // the one to determine the size of the index space
        return NumberCache(
          *std::max_element(new_numbers.begin(), new_numbers.end()) + 1);
      }



      template <class DoFHandlerType>
      NumberCache
      Sequential<DoFHandlerType>::renumber_mg_dofs(
        const unsigned int                          level,
        const std::vector<types::global_dof_index> &new_numbers) const
      {
        Implementation::renumber_mg_dofs(
          new_numbers, IndexSet(0), *dof_handler, level, true);

        // return a sequential, complete index set
        return NumberCache(new_numbers.size());
      }


      /* --------------------- class ParallelShared ---------------- */


      template <class DoFHandlerType>
      ParallelShared<DoFHandlerType>::ParallelShared(
        DoFHandlerType &dof_handler)
        : dof_handler(&dof_handler)
      {}



      namespace
      {
        template <class DoFHandlerType>
        std::vector<types::subdomain_id>
        get_dof_subdomain_association(const DoFHandlerType &        dof_handler,
                                      const types::global_dof_index n_dofs,
                                      const unsigned int            n_procs)
        {
          (void)n_procs;
          std::vector<types::subdomain_id> subdomain_association(
            n_dofs, numbers::invalid_subdomain_id);
          std::vector<types::global_dof_index> local_dof_indices;
          local_dof_indices.reserve(DoFTools::max_dofs_per_cell(dof_handler));

          // loop over all cells and record which subdomain a DoF belongs to.
          // give to the smaller subdomain_id in case it is on an interface
          typename DoFHandlerType::active_cell_iterator
            cell = dof_handler.begin_active(),
            endc = dof_handler.end();
          for (; cell != endc; ++cell)
            {
              // get the owner of the cell; note that we have made sure above
              // that all cells are either locally owned or ghosts (not
              // artificial), so this call will always yield the true owner
              const types::subdomain_id subdomain_id = cell->subdomain_id();
              const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
              local_dof_indices.resize(dofs_per_cell);
              cell->get_dof_indices(local_dof_indices);

              // set subdomain ids. if dofs already have their values set then
              // they must be on partition interfaces. In that case assign them
              // to the processor with the smaller subdomain id.
              for (unsigned int i = 0; i < dofs_per_cell; ++i)
                if (subdomain_association[local_dof_indices[i]] ==
                    numbers::invalid_subdomain_id)
                  subdomain_association[local_dof_indices[i]] = subdomain_id;
                else if (subdomain_association[local_dof_indices[i]] >
                         subdomain_id)
                  {
                    subdomain_association[local_dof_indices[i]] = subdomain_id;
                  }
            }

          Assert(std::find(subdomain_association.begin(),
                           subdomain_association.end(),
                           numbers::invalid_subdomain_id) ==
                   subdomain_association.end(),
                 ExcInternalError());

          Assert(*std::max_element(subdomain_association.begin(),
                                   subdomain_association.end()) < n_procs,
                 ExcInternalError());

          return subdomain_association;
        }


        template <class DoFHandlerType>
        std::vector<types::subdomain_id>
        get_dof_level_subdomain_association(
          const DoFHandlerType &        dof_handler,
          const types::global_dof_index n_dofs_on_level,
          const unsigned int            n_procs,
          const unsigned int            level)
        {
          (void)n_procs;
          std::vector<types::subdomain_id> level_subdomain_association(
            n_dofs_on_level, numbers::invalid_subdomain_id);
          std::vector<types::global_dof_index> local_dof_indices;
          local_dof_indices.reserve(DoFTools::max_dofs_per_cell(dof_handler));

          // loop over all cells and record which subdomain a DoF belongs to.
          // interface goes to proccessor with smaller subdomain id
          typename DoFHandlerType::cell_iterator cell =
                                                   dof_handler.begin(level),
                                                 endc = dof_handler.end(level);
          for (; cell != endc; ++cell)
            {
              // get the owner of the cell; note that we have made sure above
              // that all cells are either locally owned or ghosts (not
              // artificial), so this call will always yield the true owner
              const types::subdomain_id level_subdomain_id =
                cell->level_subdomain_id();
              const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
              local_dof_indices.resize(dofs_per_cell);
              cell->get_mg_dof_indices(local_dof_indices);

              // set level subdomain ids. if dofs already have their values set
              // then they must be on partition interfaces. In that case assign
              // them to the processor with the smaller subdomain id.
              for (unsigned int i = 0; i < dofs_per_cell; ++i)
                if (level_subdomain_association[local_dof_indices[i]] ==
                    numbers::invalid_subdomain_id)
                  level_subdomain_association[local_dof_indices[i]] =
                    level_subdomain_id;
                else if (level_subdomain_association[local_dof_indices[i]] >
                         level_subdomain_id)
                  {
                    level_subdomain_association[local_dof_indices[i]] =
                      level_subdomain_id;
                  }
            }

          Assert(std::find(level_subdomain_association.begin(),
                           level_subdomain_association.end(),
                           numbers::invalid_subdomain_id) ==
                   level_subdomain_association.end(),
                 ExcInternalError());

          Assert(*std::max_element(level_subdomain_association.begin(),
                                   level_subdomain_association.end()) < n_procs,
                 ExcInternalError());

          return level_subdomain_association;
        }
      } // namespace



      template <class DoFHandlerType>
      NumberCache
      ParallelShared<DoFHandlerType>::distribute_dofs() const
      {
        const unsigned int dim      = DoFHandlerType::dimension;
        const unsigned int spacedim = DoFHandlerType::space_dimension;

        const parallel::shared::Triangulation<dim, spacedim> *tr =
          (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
            &this->dof_handler->get_triangulation()));
        Assert(tr != nullptr, ExcInternalError());

        const unsigned int n_procs =
          Utilities::MPI::n_mpi_processes(tr->get_communicator());

        // If the underlying shared::Tria allows artificial cells,
        // then save the current set of subdomain ids, and set
        // subdomain ids to the "true" owner of each cell. we later
        // restore these flags
        std::vector<types::subdomain_id> saved_subdomain_ids;
        if (tr->with_artificial_cells())
          {
            saved_subdomain_ids.resize(tr->n_active_cells());

            typename parallel::shared::Triangulation<dim, spacedim>::
              active_cell_iterator
                cell = this->dof_handler->get_triangulation().begin_active(),
                endc = this->dof_handler->get_triangulation().end();

            const std::vector<types::subdomain_id> &true_subdomain_ids =
              tr->get_true_subdomain_ids_of_cells();

            for (unsigned int index = 0; cell != endc; ++cell, ++index)
              {
                saved_subdomain_ids[index] = cell->subdomain_id();
                cell->set_subdomain_id(true_subdomain_ids[index]);
              }
          }

        // first let the sequential algorithm do its magic. it is going to
        // enumerate DoFs on all cells, regardless of owner
        const types::global_dof_index n_dofs =
          Implementation::distribute_dofs(numbers::invalid_subdomain_id,
                                          *this->dof_handler);

        // then re-enumerate them based on their subdomain association.
        // for this, we first have to identify for each current DoF
        // index which subdomain they belong to. ideally, we would
        // like to call DoFRenumbering::subdomain_wise(), but
        // because the NumberCache of the current DoFHandler is not
        // fully set up yet, we can't quite do that. also, that
        // function has to deal with other kinds of triangulations as
        // well, whereas we here know what kind of triangulation
        // we have and can simplify the code accordingly
        std::vector<types::global_dof_index> new_dof_indices(
          n_dofs, numbers::invalid_dof_index);
        {
          // first get the association of each dof with a subdomain and
          // determine the total number of subdomain ids used
          const std::vector<types::subdomain_id> subdomain_association =
            get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);

          // then renumber the subdomains by first looking at those belonging
          // to subdomain 0, then those of subdomain 1, etc. note that the
          // algorithm is stable, i.e. if two dofs i,j have i<j and belong to
          // the same subdomain, then they will be in this order also after
          // reordering
          types::global_dof_index next_free_index = 0;
          for (types::subdomain_id subdomain = 0; subdomain < n_procs;
               ++subdomain)
            for (types::global_dof_index i = 0; i < n_dofs; ++i)
              if (subdomain_association[i] == subdomain)
                {
                  Assert(new_dof_indices[i] == numbers::invalid_dof_index,
                         ExcInternalError());
                  new_dof_indices[i] = next_free_index;
                  ++next_free_index;
                }

          // we should have numbered all dofs
          Assert(next_free_index == n_dofs, ExcInternalError());
          Assert(std::find(new_dof_indices.begin(),
                           new_dof_indices.end(),
                           numbers::invalid_dof_index) == new_dof_indices.end(),
                 ExcInternalError());
        }
        // finally do the renumbering. we can use the sequential
        // version of the function because we do things on all
        // cells and all cells have their subdomain ids and DoFs
        // correctly set
        Implementation::renumber_dofs(new_dof_indices,
                                      IndexSet(0),
                                      *this->dof_handler,
                                      true);

        // update the number cache. for this, we first have to find the
        // subdomain association for each DoF again following renumbering, from
        // which we can then compute the IndexSets of locally owned DoFs for all
        // processors. all other fields then follow from this
        //
        // given the way we enumerate degrees of freedom, the locally owned
        // ranges must all be contiguous and consecutive. this makes filling
        // the IndexSets cheap. an assertion at the top verifies that this
        // assumption is true
        const std::vector<types::subdomain_id> subdomain_association =
          get_dof_subdomain_association(*this->dof_handler, n_dofs, n_procs);

        for (unsigned int i = 1; i < n_dofs; ++i)
          Assert(subdomain_association[i] >= subdomain_association[i - 1],
                 ExcInternalError());

        std::vector<IndexSet> locally_owned_dofs_per_processor(
          n_procs, IndexSet(n_dofs));
        {
          // we know that the set of subdomain indices is contiguous from
          // the assertion above; find the start and end index for each
          // processor, taking into account that sometimes a processor
          // may not in fact have any DoFs at all. we do the latter
          // by just identifying contiguous ranges of subdomain_ids
          // and filling IndexSets for those subdomains; subdomains
          // that don't appear will lead to IndexSets that are simply
          // never touched and remain empty as initialized above.
          unsigned int start_index = 0;
          unsigned int end_index   = 0;
          while (start_index < n_dofs)
            {
              while ((end_index) < n_dofs &&
                     (subdomain_association[end_index] ==
                      subdomain_association[start_index]))
                ++end_index;

              // we've now identified a range of same indices. set that
              // range in the corresponding IndexSet
              if (end_index > start_index)
                {
                  const unsigned int subdomain_owner =
                    subdomain_association[start_index];
                  locally_owned_dofs_per_processor[subdomain_owner].add_range(
                    start_index, end_index);
                }

              // then move on to thinking about the next range
              start_index = end_index;
            }
        }

        // finally, restore current subdomain ids
        if (tr->with_artificial_cells())
          {
            typename parallel::shared::Triangulation<dim, spacedim>::
              active_cell_iterator
                cell = this->dof_handler->get_triangulation().begin_active(),
                endc = this->dof_handler->get_triangulation().end();

            for (unsigned int index = 0; cell != endc; ++cell, ++index)
              cell->set_subdomain_id(saved_subdomain_ids[index]);
          }

        // return a NumberCache object made up from the sets of locally
        // owned DoFs
        return NumberCache(
          locally_owned_dofs_per_processor,
          this->dof_handler->get_triangulation().locally_owned_subdomain());
      }



      template <class DoFHandlerType>
      std::vector<NumberCache>
      ParallelShared<DoFHandlerType>::distribute_mg_dofs() const
      {
        const unsigned int dim      = DoFHandlerType::dimension;
        const unsigned int spacedim = DoFHandlerType::space_dimension;

        const parallel::shared::Triangulation<dim, spacedim> *tr =
          (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
            &this->dof_handler->get_triangulation()));
        Assert(tr != nullptr, ExcInternalError());

        const unsigned int n_procs =
          Utilities::MPI::n_mpi_processes(tr->get_communicator());
        const unsigned int n_levels = tr->n_global_levels();

        std::vector<NumberCache> number_caches;
        number_caches.reserve(n_levels);

        // We create an index set for each level
        for (unsigned int lvl = 0; lvl < n_levels; ++lvl)
          {
            // If the underlying shared::Tria allows artificial cells,
            // then save the current set of level subdomain ids, and set
            // subdomain ids to the "true" owner of each cell. we later
            // restore these flags
            // Note: "allows_artificial_cells" is currently enforced for
            //       MG computations.
            std::vector<types::subdomain_id> saved_level_subdomain_ids;
            saved_level_subdomain_ids.resize(tr->n_cells(lvl));
            {
              typename parallel::shared::Triangulation<dim,
                                                       spacedim>::cell_iterator
                cell = this->dof_handler->get_triangulation().begin(lvl),
                endc = this->dof_handler->get_triangulation().end(lvl);

              const std::vector<types::subdomain_id> &true_level_subdomain_ids =
                tr->get_true_level_subdomain_ids_of_cells(lvl);

              for (unsigned int index = 0; cell != endc; ++cell, ++index)
                {
                  saved_level_subdomain_ids[index] = cell->level_subdomain_id();
                  cell->set_level_subdomain_id(true_level_subdomain_ids[index]);
                }
            }

            // Next let the sequential algorithm do its magic. it is going to
            // enumerate DoFs on all cells on the given level, regardless of
            // owner
            const types::global_dof_index n_dofs_on_level =
              Implementation::distribute_dofs_on_level(
                numbers::invalid_subdomain_id, *this->dof_handler, lvl);

            // then re-enumerate them based on their level subdomain
            // association. for this, we first have to identify for each current
            // DoF index which subdomain they belong to. ideally, we would like
            // to call DoFRenumbering::subdomain_wise(), but because the
            // NumberCache of the current DoFHandler is not fully set up yet, we
            // can't quite do that. also, that function has to deal with other
            // kinds of triangulations as well, whereas we here know what kind
            // of triangulation we have and can simplify the code accordingly
            std::vector<types::global_dof_index> new_dof_indices(
              n_dofs_on_level, numbers::invalid_dof_index);
            {
              // first get the association of each dof with a subdomain and
              // determine the total number of subdomain ids used
              const std::vector<types::subdomain_id>
                level_subdomain_association =
                  get_dof_level_subdomain_association(*this->dof_handler,
                                                      n_dofs_on_level,
                                                      n_procs,
                                                      lvl);

              // then renumber the subdomains by first looking at those
              // belonging to subdomain 0, then those of subdomain 1, etc. note
              // that the algorithm is stable, i.e. if two dofs i,j have i<j and
              // belong to the same subdomain, then they will be in this order
              // also after reordering
              types::global_dof_index next_free_index = 0;
              for (types::subdomain_id level_subdomain = 0;
                   level_subdomain < n_procs;
                   ++level_subdomain)
                for (types::global_dof_index i = 0; i < n_dofs_on_level; ++i)
                  if (level_subdomain_association[i] == level_subdomain)
                    {
                      Assert(new_dof_indices[i] == numbers::invalid_dof_index,
                             ExcInternalError());
                      new_dof_indices[i] = next_free_index;
                      ++next_free_index;
                    }

              // we should have numbered all dofs
              Assert(next_free_index == n_dofs_on_level, ExcInternalError());
              Assert(std::find(new_dof_indices.begin(),
                               new_dof_indices.end(),
                               numbers::invalid_dof_index) ==
                       new_dof_indices.end(),
                     ExcInternalError());
            }

            // finally do the renumbering. we can use the sequential
            // version of the function because we do things on all
            // cells and all cells have their subdomain ids and DoFs
            // correctly set
            Implementation::renumber_mg_dofs(
              new_dof_indices, IndexSet(0), *this->dof_handler, lvl, true);

            // update the number cache. for this, we first have to find the
            // level subdomain association for each DoF again following
            // renumbering, from which we can then compute the IndexSets of
            // locally owned DoFs for all processors. all other fields then
            // follow from this
            //
            // given the way we enumerate degrees of freedom, the locally owned
            // ranges must all be contiguous and consecutive. this makes filling
            // the IndexSets cheap. an assertion at the top verifies that this
            // assumption is true
            const std::vector<types::subdomain_id> level_subdomain_association =
              get_dof_level_subdomain_association(*this->dof_handler,
                                                  n_dofs_on_level,
                                                  n_procs,
                                                  lvl);

            for (unsigned int i = 1; i < n_dofs_on_level; ++i)
              Assert(level_subdomain_association[i] >=
                       level_subdomain_association[i - 1],
                     ExcInternalError());

            std::vector<IndexSet> locally_owned_dofs_per_processor(
              n_procs, IndexSet(n_dofs_on_level));
            {
              // we know that the set of subdomain indices is contiguous from
              // the assertion above; find the start and end index for each
              // processor, taking into account that sometimes a processor
              // may not in fact have any DoFs at all. we do the latter
              // by just identifying contiguous ranges of level_subdomain_ids
              // and filling IndexSets for those subdomains; subdomains
              // that don't appear will lead to IndexSets that are simply
              // never touched and remain empty as initialized above.
              unsigned int start_index = 0;
              unsigned int end_index   = 0;
              while (start_index < n_dofs_on_level)
                {
                  while ((end_index) < n_dofs_on_level &&
                         (level_subdomain_association[end_index] ==
                          level_subdomain_association[start_index]))
                    ++end_index;

                  // we've now identified a range of same indices. set that
                  // range in the corresponding IndexSet
                  if (end_index > start_index)
                    {
                      const unsigned int level_subdomain_owner =
                        level_subdomain_association[start_index];
                      locally_owned_dofs_per_processor[level_subdomain_owner]
                        .add_range(start_index, end_index);
                    }

                  // then move on to thinking about the next range
                  start_index = end_index;
                }
            }

            // finally, restore current level subdomain ids
            {
              typename parallel::shared::Triangulation<dim,
                                                       spacedim>::cell_iterator
                cell = this->dof_handler->get_triangulation().begin(lvl),
                endc = this->dof_handler->get_triangulation().end(lvl);

              for (unsigned int index = 0; cell != endc; ++cell, ++index)
                cell->set_level_subdomain_id(saved_level_subdomain_ids[index]);

              // add NumberCache for current level
              number_caches.emplace_back(
                NumberCache(locally_owned_dofs_per_processor,
                            this->dof_handler->get_triangulation()
                              .locally_owned_subdomain()));
            }
          }

        return number_caches;
      }



      template <class DoFHandlerType>
      NumberCache
      ParallelShared<DoFHandlerType>::renumber_dofs(
        const std::vector<types::global_dof_index> &new_numbers) const
      {
#ifndef DEAL_II_WITH_MPI
        (void)new_numbers;
        Assert(false, ExcNotImplemented());
        return NumberCache();
#else
        const unsigned int dim      = DoFHandlerType::dimension;
        const unsigned int spacedim = DoFHandlerType::space_dimension;

        // Similar to distribute_dofs() we need to have a special treatment in
        // case artificial cells are present.
        const parallel::shared::Triangulation<dim, spacedim> *tr =
          (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
            &this->dof_handler->get_triangulation()));
        Assert(tr != nullptr, ExcInternalError());

        typename parallel::shared::Triangulation<dim,
                                                 spacedim>::active_cell_iterator
          cell = this->dof_handler->get_triangulation().begin_active(),
          endc = this->dof_handler->get_triangulation().end();
        std::vector<types::subdomain_id> current_subdomain_ids(
          tr->n_active_cells());
        const std::vector<types::subdomain_id> &true_subdomain_ids =
          tr->get_true_subdomain_ids_of_cells();
        if (tr->with_artificial_cells())
          for (unsigned int index = 0; cell != endc; cell++, index++)
            {
              current_subdomain_ids[index] = cell->subdomain_id();
              cell->set_subdomain_id(true_subdomain_ids[index]);
            }

        std::vector<types::global_dof_index> global_gathered_numbers(
          this->dof_handler->n_dofs(), 0);
        // as we call DoFRenumbering::subdomain_wise (*dof_handler) from
        // distribute_dofs(), we need to support sequential-like input.
        // Distributed-like input from, for example, component_wise renumbering
        // is also supported.
        if (new_numbers.size() == this->dof_handler->n_dofs())
          {
            global_gathered_numbers = new_numbers;
          }
        else
          {
            Assert(new_numbers.size() ==
                     this->dof_handler->locally_owned_dofs().n_elements(),
                   ExcInternalError());
            const unsigned int n_cpu =
              Utilities::MPI::n_mpi_processes(tr->get_communicator());
            std::vector<types::global_dof_index> gathered_new_numbers(
              this->dof_handler->n_dofs(), 0);
            Assert(Utilities::MPI::this_mpi_process(tr->get_communicator()) ==
                     this->dof_handler->get_triangulation()
                       .locally_owned_subdomain(),
                   ExcInternalError())

            // gather new numbers among processors into one vector
            {
              std::vector<types::global_dof_index> new_numbers_copy(
                new_numbers);

              // store the number of elements that are to be received from each
              // process
              std::vector<int> rcounts(n_cpu);

              types::global_dof_index shift = 0;
              // set rcounts based on new_numbers:
              int cur_count = new_numbers_copy.size();
              int ierr      = MPI_Allgather(&cur_count,
                                       1,
                                       MPI_INT,
                                       rcounts.data(),
                                       1,
                                       MPI_INT,
                                       tr->get_communicator());
              AssertThrowMPI(ierr);

              // compute the displacements (relative to recvbuf)
              // at which to place the incoming data from process i
              std::vector<int> displacements(n_cpu);
              for (unsigned int i = 0; i < n_cpu; i++)
                {
                  displacements[i] = shift;
                  shift += rcounts[i];
                }
              Assert(((int)new_numbers_copy.size()) ==
                       rcounts[Utilities::MPI::this_mpi_process(
                         tr->get_communicator())],
                     ExcInternalError());
              ierr = MPI_Allgatherv(new_numbers_copy.data(),
                                    new_numbers_copy.size(),
                                    DEAL_II_DOF_INDEX_MPI_TYPE,
                                    gathered_new_numbers.data(),
                                    rcounts.data(),
                                    displacements.data(),
                                    DEAL_II_DOF_INDEX_MPI_TYPE,
                                    tr->get_communicator());
              AssertThrowMPI(ierr);
            }

            // put new numbers according to the current
            // locally_owned_dofs_per_processor IndexSets
            types::global_dof_index shift = 0;
            // flag_1 and flag_2 are
            // used to control that there is a
            // one-to-one relation between old and new DoFs.
            std::vector<unsigned int> flag_1(this->dof_handler->n_dofs(), 0);
            std::vector<unsigned int> flag_2(this->dof_handler->n_dofs(), 0);
            for (unsigned int i = 0; i < n_cpu; i++)
              {
                const IndexSet iset =
                  this->dof_handler->locally_owned_dofs_per_processor()[i];
                for (types::global_dof_index ind = 0; ind < iset.n_elements();
                     ind++)
                  {
                    const types::global_dof_index target =
                      iset.nth_index_in_set(ind);
                    const types::global_dof_index value =
                      gathered_new_numbers[shift + ind];
                    Assert(target < this->dof_handler->n_dofs(),
                           ExcInternalError());
                    Assert(value < this->dof_handler->n_dofs(),
                           ExcInternalError());
                    global_gathered_numbers[target] = value;
                    flag_1[target]++;
                    flag_2[value]++;
                  }
                shift += iset.n_elements();
              }

            Assert(*std::max_element(flag_1.begin(), flag_1.end()) == 1,
                   ExcInternalError());
            Assert(*std::min_element(flag_1.begin(), flag_1.end()) == 1,
                   ExcInternalError());
            Assert((*std::max_element(flag_2.begin(), flag_2.end())) == 1,
                   ExcInternalError());
            Assert((*std::min_element(flag_2.begin(), flag_2.end())) == 1,
                   ExcInternalError());
          }

        // let the sequential algorithm do its magic; ignore the
        // return type, but reconstruct the number cache based on
        // which DoFs each process owns
        Implementation::renumber_dofs(global_gathered_numbers,
                                      IndexSet(0),
                                      *this->dof_handler,
                                      true);

        const NumberCache number_cache(
          DoFTools::locally_owned_dofs_per_subdomain(*this->dof_handler),
          this->dof_handler->get_triangulation().locally_owned_subdomain());

        // restore artificial cells
        cell = tr->begin_active();
        if (tr->with_artificial_cells())
          for (unsigned int index = 0; cell != endc; cell++, index++)
            cell->set_subdomain_id(current_subdomain_ids[index]);

        return number_cache;
#endif
      }



      template <class DoFHandlerType>
      NumberCache
      ParallelShared<DoFHandlerType>::renumber_mg_dofs(
        const unsigned int /*level*/,
        const std::vector<types::global_dof_index> & /*new_numbers*/) const
      {
        // multigrid is not currently implemented for shared triangulations
        Assert(false, ExcNotImplemented());

        return NumberCache();
      }



      /* --------------------- class ParallelDistributed ---------------- */

#ifdef DEAL_II_WITH_P4EST

      namespace
      {
        template <int dim>
        struct CellDataTransferBuffer
        {
          std::vector<unsigned int> tree_indices;
          std::vector<typename ::internal::p4est::types<dim>::quadrant>
                                                       quadrants;
          std::vector<::types::global_dof_index> dof_numbers_and_indices;


          template <class Archive>
          void
          save(Archive &ar, const unsigned int /*version*/) const
          {
            // we would like to directly serialize the 'quadrants' vector,
            // but the element type is internal to p4est and does not
            // know how to serialize itself. consequently, first copy it over
            // to an array of bytes, and then serialize that
            std::vector<char> quadrants_as_chars(sizeof(quadrants[0]) *
                                                 quadrants.size());
            if (quadrants_as_chars.size() > 0)
              {
                Assert(quadrants.data() != nullptr, ExcInternalError());
                std::memcpy(quadrants_as_chars.data(),
                            quadrants.data(),
                            quadrants_as_chars.size());
              }

            // now serialize everything
            ar &quadrants_as_chars &tree_indices &dof_numbers_and_indices;
          }

          template <class Archive>
          void
          load(Archive &ar, const unsigned int /*version*/)
          {
            // undo the copying trick from the 'save' function
            std::vector<char> quadrants_as_chars;
            ar &quadrants_as_chars &tree_indices &dof_numbers_and_indices;

            if (quadrants_as_chars.size() > 0)
              {
                quadrants.resize(quadrants_as_chars.size() /
                                 sizeof(quadrants[0]));
                std::memcpy(quadrants.data(),
                            quadrants_as_chars.data(),
                            quadrants_as_chars.size());
              }
            else
              quadrants.clear();
          }

          BOOST_SERIALIZATION_SPLIT_MEMBER()


          
          std::vector<char>
          pack_data() const
          {
            // set up a buffer and then use it as the target of a compressing
            // stream into which we serialize the current object
            std::vector<char> buffer;
            {
#  ifdef DEAL_II_WITH_ZLIB
              boost::iostreams::filtering_ostream out;
              out.push(
                boost::iostreams::gzip_compressor(boost::iostreams::gzip_params(
                  boost::iostreams::gzip::best_compression)));
              out.push(boost::iostreams::back_inserter(buffer));

              boost::archive::binary_oarchive archive(out);

              archive << *this;
              out.flush();
#  else
              std::ostringstream              out;
              boost::archive::binary_oarchive archive(out);
              archive << *this;
              const std::string &s = out.str();
              buffer.reserve(s.size());
              buffer.assign(s.begin(), s.end());
#  endif
            }

            return buffer;
          }


          void
          unpack_data(const std::vector<char> &buffer)
          {
            std::string decompressed_buffer;

            // first decompress the buffer
            {
#  ifdef DEAL_II_WITH_ZLIB
              boost::iostreams::filtering_ostream decompressing_stream;
              decompressing_stream.push(boost::iostreams::gzip_decompressor());
              decompressing_stream.push(
                boost::iostreams::back_inserter(decompressed_buffer));
              decompressing_stream.write(buffer.data(), buffer.size());
#  else
              decompressed_buffer.assign(buffer.begin(), buffer.end());
#  endif
            }

            // then restore the object from the buffer
            std::istringstream              in(decompressed_buffer);
            boost::archive::binary_iarchive archive(in);

            archive >> *this;
          }
        };



        template <int dim, int spacedim>
        void
        get_mg_dofindices_recursively(
          const parallel::distributed::Triangulation<dim, spacedim> &tria,
          const typename ::internal::p4est::types<dim>::quadrant
            &p4est_cell,
          const typename DoFHandler<dim, spacedim>::level_cell_iterator
            &dealii_cell,
          const typename ::internal::p4est::types<dim>::quadrant
            &                          quadrant,
          CellDataTransferBuffer<dim> &cell_data_transfer_buffer)
        {
          if (internal::p4est::quadrant_is_equal<dim>(p4est_cell, quadrant))
            {
              // why would somebody request a cell that is not ours?
              Assert(dealii_cell->level_subdomain_id() ==
                       tria.locally_owned_subdomain(),
                     ExcInternalError());


              std::vector<::types::global_dof_index> local_dof_indices(
                dealii_cell->get_fe().dofs_per_cell);
              dealii_cell->get_mg_dof_indices(local_dof_indices);

              cell_data_transfer_buffer.dof_numbers_and_indices.push_back(
                dealii_cell->get_fe().dofs_per_cell);
              cell_data_transfer_buffer.dof_numbers_and_indices.insert(
                cell_data_transfer_buffer.dof_numbers_and_indices.end(),
                local_dof_indices.begin(),
                local_dof_indices.end());
              return; // we are done
            }

          if (!dealii_cell->has_children())
            return;

          if (!internal::p4est::quadrant_is_ancestor<dim>(p4est_cell, quadrant))
            return;

          typename ::internal::p4est::types<dim>::quadrant
            p4est_child[GeometryInfo<dim>::max_children_per_cell];
          internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);

          for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
               ++c)
            get_mg_dofindices_recursively<dim, spacedim>(
              tria,
              p4est_child[c],
              dealii_cell->child(c),
              quadrant,
              cell_data_transfer_buffer);
        }


        template <int dim, int spacedim>
        void
        find_marked_mg_ghost_cells_recursively(
          const typename parallel::distributed::Triangulation<dim, spacedim>
            &                tria,
          const unsigned int tree_index,
          const typename DoFHandler<dim, spacedim>::level_cell_iterator
            &dealii_cell,
          const typename ::internal::p4est::types<dim>::quadrant
            &p4est_cell,
          std::map<::types::subdomain_id, CellDataTransferBuffer<dim>>
            &neighbor_cell_list)
        {
          // recurse...
          if (dealii_cell->has_children())
            {
              typename ::internal::p4est::types<dim>::quadrant
                p4est_child[GeometryInfo<dim>::max_children_per_cell];
              internal::p4est::init_quadrant_children<dim>(p4est_cell,
                                                           p4est_child);


              for (unsigned int c = 0;
                   c < GeometryInfo<dim>::max_children_per_cell;
                   ++c)
                find_marked_mg_ghost_cells_recursively<dim, spacedim>(
                  tria,
                  tree_index,
                  dealii_cell->child(c),
                  p4est_child[c],
                  neighbor_cell_list);
            }

          if (dealii_cell->user_flag_set() &&
              dealii_cell->level_subdomain_id() !=
                tria.locally_owned_subdomain())
            {
              neighbor_cell_list[dealii_cell->level_subdomain_id()]
                .tree_indices.push_back(tree_index);
              neighbor_cell_list[dealii_cell->level_subdomain_id()]
                .quadrants.push_back(p4est_cell);
            }
        }


        template <int dim, int spacedim>
        void
        set_mg_dofindices_recursively(
          const parallel::distributed::Triangulation<dim, spacedim> &tria,
          const typename ::internal::p4est::types<dim>::quadrant
            &p4est_cell,
          const typename DoFHandler<dim, spacedim>::level_cell_iterator
            &dealii_cell,
          const typename ::internal::p4est::types<dim>::quadrant
            &                              quadrant,
          ::types::global_dof_index *dofs)
        {
          if (internal::p4est::quadrant_is_equal<dim>(p4est_cell, quadrant))
            {
              Assert(dealii_cell->level_subdomain_id() !=
                       ::numbers::artificial_subdomain_id,
                     ExcInternalError());

              // update dof indices of cell
              std::vector<::types::global_dof_index> dof_indices(
                dealii_cell->get_fe().dofs_per_cell);
              dealii_cell->get_mg_dof_indices(dof_indices);

              bool complete = true;
              for (unsigned int i = 0; i < dof_indices.size(); ++i)
                if (dofs[i] != numbers::invalid_dof_index)
                  {
                    Assert((dof_indices[i] == (numbers::invalid_dof_index)) ||
                             (dof_indices[i] == dofs[i]),
                           ExcInternalError());
                    dof_indices[i] = dofs[i];
                  }
                else
                  complete = false;

              if (!complete)
                const_cast<
                  typename DoFHandler<dim, spacedim>::level_cell_iterator &>(
                  dealii_cell)
                  ->set_user_flag();
              else
                const_cast<
                  typename DoFHandler<dim, spacedim>::level_cell_iterator &>(
                  dealii_cell)
                  ->clear_user_flag();

              const_cast<
                typename DoFHandler<dim, spacedim>::level_cell_iterator &>(
                dealii_cell)
                ->set_mg_dof_indices(dof_indices);
              return;
            }

          if (!dealii_cell->has_children())
            return;

          if (!internal::p4est::quadrant_is_ancestor<dim>(p4est_cell, quadrant))
            return;

          typename ::internal::p4est::types<dim>::quadrant
            p4est_child[GeometryInfo<dim>::max_children_per_cell];
          internal::p4est::init_quadrant_children<dim>(p4est_cell, p4est_child);

          for (unsigned int c = 0; c < GeometryInfo<dim>::max_children_per_cell;
               ++c)
            set_mg_dofindices_recursively<dim, spacedim>(
              tria, p4est_child[c], dealii_cell->child(c), quadrant, dofs);
        }



        template <int dim, int spacedim, class DoFHandlerType>
        void
        communicate_mg_ghost_cells(
          const typename parallel::distributed::Triangulation<dim, spacedim>
            &             tria,
          DoFHandlerType &dof_handler,
          const std::vector<::types::global_dof_index>
            &coarse_cell_to_p4est_tree_permutation,
          const std::vector<::types::global_dof_index>
            &p4est_tree_to_coarse_cell_permutation)
        {
          // build list of cells to request for each neighbor
          std::set<::types::subdomain_id> level_ghost_owners =
            tria.level_ghost_owners();
          using cellmap_t =
            std::map<::types::subdomain_id, CellDataTransferBuffer<dim>>;
          cellmap_t neighbor_cell_list;
          for (std::set<::types::subdomain_id>::iterator it =
                 level_ghost_owners.begin();
               it != level_ghost_owners.end();
               ++it)
            neighbor_cell_list.insert(
              std::make_pair(*it, CellDataTransferBuffer<dim>()));

          for (typename DoFHandlerType::level_cell_iterator cell =
                 dof_handler.begin(0);
               cell != dof_handler.end(0);
               ++cell)
            {
              typename ::internal::p4est::types<dim>::quadrant
                p4est_coarse_cell;
              internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);

              find_marked_mg_ghost_cells_recursively<dim, spacedim>(
                tria,
                coarse_cell_to_p4est_tree_permutation[cell->index()],
                cell,
                p4est_coarse_cell,
                neighbor_cell_list);
            }
          Assert(level_ghost_owners.size() == neighbor_cell_list.size(),
                 ExcInternalError());

          //* send our requests:
          std::vector<std::vector<char>> sendbuffers(level_ghost_owners.size());
          std::vector<MPI_Request>       requests(level_ghost_owners.size());

          unsigned int idx = 0;
          for (typename cellmap_t::iterator it = neighbor_cell_list.begin();
               it != neighbor_cell_list.end();
               ++it, ++idx)
            {
              // pack all the data into the buffer for this recipient
              // and send it. keep data around till we can make sure
              // that the packet has been received
              sendbuffers[idx] = it->second.pack_data();
              const int ierr   = MPI_Isend(sendbuffers[idx].data(),
                                         sendbuffers[idx].size(),
                                         MPI_BYTE,
                                         it->first,
                                         10101,
                                         tria.get_communicator(),
                                         &requests[idx]);
              AssertThrowMPI(ierr);
            }

          //* receive requests and reply
          std::vector<std::vector<char>> reply_buffers(
            level_ghost_owners.size());
          std::vector<MPI_Request> reply_requests(level_ghost_owners.size());

          for (unsigned int idx = 0; idx < level_ghost_owners.size(); ++idx)
            {
              std::vector<char>           receive;
              CellDataTransferBuffer<dim> cell_data_transfer_buffer;

              MPI_Status status;
              int        len;
              int        ierr = MPI_Probe(MPI_ANY_SOURCE,
                                   10101,
                                   tria.get_communicator(),
                                   &status);
              AssertThrowMPI(ierr);
              ierr = MPI_Get_count(&status, MPI_BYTE, &len);
              AssertThrowMPI(ierr);
              receive.resize(len);

              char *ptr = receive.data();
              ierr      = MPI_Recv(ptr,
                              len,
                              MPI_BYTE,
                              status.MPI_SOURCE,
                              status.MPI_TAG,
                              tria.get_communicator(),
                              &status);
              AssertThrowMPI(ierr);

              cell_data_transfer_buffer.unpack_data(receive);

              // store the dof indices for each cell
              for (unsigned int c = 0;
                   c < cell_data_transfer_buffer.tree_indices.size();
                   ++c)
                {
                  typename DoFHandlerType::level_cell_iterator cell(
                    &dof_handler.get_triangulation(),
                    0,
                    p4est_tree_to_coarse_cell_permutation
                      [cell_data_transfer_buffer.tree_indices[c]],
                    &dof_handler);

                  typename ::internal::p4est::types<dim>::quadrant
                    p4est_coarse_cell;
                  internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);

                  get_mg_dofindices_recursively<dim, spacedim>(
                    tria,
                    p4est_coarse_cell,
                    cell,
                    cell_data_transfer_buffer.quadrants[c],
                    cell_data_transfer_buffer);
                }

              // send reply
              reply_buffers[idx] = cell_data_transfer_buffer.pack_data();
              ierr               = MPI_Isend(&(reply_buffers[idx])[0],
                               reply_buffers[idx].size(),
                               MPI_BYTE,
                               status.MPI_SOURCE,
                               10102,
                               tria.get_communicator(),
                               &reply_requests[idx]);
              AssertThrowMPI(ierr);
            }

          //* finally receive the replies
          for (unsigned int idx = 0; idx < level_ghost_owners.size(); ++idx)
            {
              std::vector<char>           receive;
              CellDataTransferBuffer<dim> cell_data_transfer_buffer;

              MPI_Status status;
              int        len;
              int        ierr = MPI_Probe(MPI_ANY_SOURCE,
                                   10102,
                                   tria.get_communicator(),
                                   &status);
              AssertThrowMPI(ierr);
              ierr = MPI_Get_count(&status, MPI_BYTE, &len);
              AssertThrowMPI(ierr);
              receive.resize(len);

              char *ptr = receive.data();
              ierr      = MPI_Recv(ptr,
                              len,
                              MPI_BYTE,
                              status.MPI_SOURCE,
                              status.MPI_TAG,
                              tria.get_communicator(),
                              &status);
              AssertThrowMPI(ierr);

              cell_data_transfer_buffer.unpack_data(receive);
              if (cell_data_transfer_buffer.tree_indices.size() == 0)
                continue;

              // set the dof indices for each cell
              ::types::global_dof_index *dofs =
                cell_data_transfer_buffer.dof_numbers_and_indices.data();
              for (unsigned int c = 0;
                   c < cell_data_transfer_buffer.tree_indices.size();
                   ++c, dofs += 1 + dofs[0])
                {
                  typename DoFHandlerType::level_cell_iterator cell(
                    &tria,
                    0,
                    p4est_tree_to_coarse_cell_permutation
                      [cell_data_transfer_buffer.tree_indices[c]],
                    &dof_handler);

                  typename ::internal::p4est::types<dim>::quadrant
                    p4est_coarse_cell;
                  internal::p4est::init_coarse_quadrant<dim>(p4est_coarse_cell);

                  Assert(cell->get_fe().dofs_per_cell == dofs[0],
                         ExcInternalError());

                  set_mg_dofindices_recursively<dim, spacedim>(
                    tria,
                    p4est_coarse_cell,
                    cell,
                    cell_data_transfer_buffer.quadrants[c],
                    dofs + 1);
                }
            }

          // complete all sends, so that we can safely destroy the
          // buffers.
          if (requests.size() > 0)
            {
              const int ierr = MPI_Waitall(requests.size(),
                                           requests.data(),
                                           MPI_STATUSES_IGNORE);
              AssertThrowMPI(ierr);
            }
          if (reply_requests.size() > 0)
            {
              const int ierr = MPI_Waitall(reply_requests.size(),
                                           reply_requests.data(),
                                           MPI_STATUSES_IGNORE);
              AssertThrowMPI(ierr);
            }
        }



        template <int spacedim>
        void
        communicate_mg_ghost_cells(
          const typename parallel::distributed::Triangulation<1, spacedim> &,
          DoFHandler<1, spacedim> &,
          const std::vector<::types::global_dof_index> &,
          const std::vector<::types::global_dof_index> &)
        {
          Assert(false, ExcNotImplemented());
        }



        template <int spacedim>
        void
        communicate_mg_ghost_cells(
          const typename parallel::distributed::Triangulation<1, spacedim> &,
          hp::DoFHandler<1, spacedim> &,
          const std::vector<::types::global_dof_index> &,
          const std::vector<::types::global_dof_index> &)
        {
          Assert(false, ExcNotImplemented());
        }



        template <int spacedim>
        void
        communicate_dof_indices_on_marked_cells(
          const DoFHandler<1, spacedim> &,
          const std::map<unsigned int, std::set<::types::subdomain_id>> &,
          const std::vector<::types::global_dof_index> &,
          const std::vector<::types::global_dof_index> &)
        {
          Assert(false, ExcNotImplemented());
        }



        template <int spacedim>
        void
        communicate_dof_indices_on_marked_cells(
          const hp::DoFHandler<1, spacedim> &,
          const std::map<unsigned int, std::set<::types::subdomain_id>> &,
          const std::vector<::types::global_dof_index> &,
          const std::vector<::types::global_dof_index> &)
        {
          Assert(false, ExcNotImplemented());
        }



        template <class DoFHandlerType>
        void
        communicate_dof_indices_on_marked_cells(
          const DoFHandlerType &dof_handler,
          const std::map<unsigned int, std::set<::types::subdomain_id>> &,
          const std::vector<::types::global_dof_index> &,
          const std::vector<::types::global_dof_index> &)
        {
#  ifndef DEAL_II_WITH_MPI
          (void)vertices_with_ghost_neighbors;
          Assert(false, ExcNotImplemented());
#  else
          const unsigned int dim = DoFHandlerType::dimension;
          const unsigned int spacedim = DoFHandlerType::space_dimension;

          // define functions that pack data on cells that are ghost cells
          // somewhere else, and unpack data on cells where we get information
          // from elsewhere
          auto pack =
            [](const typename DoFHandlerType::active_cell_iterator &cell)
            -> boost::optional<std::vector<types::global_dof_index>> {
            Assert(cell->is_locally_owned(), ExcInternalError());

            // first see whether we need to do anything at all on this cell.
            // this is determined by whether the user_flag is set on the
            // cell that indicates that the *complete* set of DoF indices
            // has not been sent
            if (cell->user_flag_set())
              {
                // get dof indices for the current cell
                std::vector<types::global_dof_index> local_dof_indices(
                  cell->get_fe().dofs_per_cell);
                cell->get_dof_indices(local_dof_indices);

                // now see if there are dof indices that were previously
                // unknown. this can only happen in phase 1, and in
                // that case we know that the user flag must have been set
                //
                // in any case, if the cell *is* complete, we do not
                // need to send the data any more in the next phase. indicate
                // this by removing the user flag
                if (std::find(local_dof_indices.begin(),
                              local_dof_indices.end(),
                              numbers::invalid_dof_index) !=
                    local_dof_indices.end())
                  {
                    Assert(cell->user_flag_set(), ExcInternalError());
                  }
                else
                  cell->clear_user_flag();

                return local_dof_indices;
              }
            else
              {
                // the fact that the user flag wasn't set means that there is
                // nothing we need to send that hasn't been sent so far.
                // so return an empty array, but also verify that indeed
                // the cell is complete
#    ifdef DEBUG
                std::vector<types::global_dof_index> local_dof_indices(
                  cell->get_fe().dofs_per_cell);
                cell->get_dof_indices(local_dof_indices);

                const bool is_complete =
                  (std::find(local_dof_indices.begin(),
                             local_dof_indices.end(),
                             numbers::invalid_dof_index) ==
                   local_dof_indices.end());
                Assert(is_complete, ExcInternalError());
#    endif
                return boost::optional<std::vector<types::global_dof_index>>();
              }
          };

          auto unpack =
            [](const typename DoFHandlerType::active_cell_iterator &cell,
               const std::vector<types::global_dof_index> &received_dof_indices)
            -> void {
            // this function should only be called on ghost cells, and
            // on top of that, only on cells that have not been
            // completed -- which we indicate via the user flag.
            // check both
            Assert(cell->is_ghost(), ExcInternalError());
            Assert(cell->user_flag_set(), ExcInternalError());

            // if we just got an incomplete array of DoF indices, then we must
            // be in the first ghost exchange and the user flag must have been
            // set. we tested that already above.
            //
            // if we did get a complete array, then we may be in the first
            // or second ghost exchange, but in any case we need not exchange
            // another time. so delete the user flag
            const bool is_complete = (std::find(received_dof_indices.begin(),
                                                received_dof_indices.end(),
                                                numbers::invalid_dof_index) ==
                                      received_dof_indices.end());
            if (is_complete)
              cell->clear_user_flag();

            // in any case, set the DoF indices on this cell. some
            // of the ones we received may still be invalid because
            // the sending processor did not know them yet, so we
            // need to merge the ones we get with those that are
            // already set here and may have already been known. for
            // those that we already know *and* get, they must obviously
            // agree
            //
            // before getting the local dof indices, we need to update the
            // cell dof indices cache because we may have set dof indices
            // on a neighboring ghost cell before this one, which may have
            // affected the dof indices we know about the current cell
            std::vector<types::global_dof_index> local_dof_indices(
              cell->get_fe().dofs_per_cell);
            cell->update_cell_dof_indices_cache();
            cell->get_dof_indices(local_dof_indices);

            for (unsigned int i = 0; i < local_dof_indices.size(); ++i)
              if (local_dof_indices[i] == numbers::invalid_dof_index)
                local_dof_indices[i] = received_dof_indices[i];
              else
                // we already know the dof index. check that there
                // is no conflict
                Assert((received_dof_indices[i] ==
                        numbers::invalid_dof_index) ||
                         (received_dof_indices[i] == local_dof_indices[i]),
                       ExcInternalError());

            const_cast<typename DoFHandlerType::active_cell_iterator &>(cell)
              ->set_dof_indices(local_dof_indices);
          };

          GridTools::exchange_cell_data_to_ghosts<
            std::vector<types::global_dof_index>,
            DoFHandlerType>(dof_handler, pack, unpack);

          // finally update the cell DoF indices caches to make sure
          // our internal data structures are consistent
          update_all_active_cell_dof_indices_caches(dof_handler);


          // have a barrier so that sends between two calls to this
          // function are not mixed up.
          //
          // this is necessary because above we just see if there are
          // messages and then receive them, without discriminating
          // where they come from and whether they were sent in phase
          // 1 or 2 (the function is called twice in a row). the need
          // for a global communication step like this barrier could
          // be avoided by receiving messages specifically from those
          // processors from which we expect messages, and by using
          // different tags for phase 1 and 2, but the cost of a
          // barrier is negligible compared to everything else we do
          // here
          if (const auto *triangulation = dynamic_cast<
                const parallel::distributed::Triangulation<dim, spacedim> *>(
                &dof_handler.get_triangulation()))
            {
              const int ierr = MPI_Barrier(triangulation->get_communicator());
              AssertThrowMPI(ierr);
            }
          else
            {
              Assert(false,
                     ExcMessage(
                       "The function communicate_dof_indices_on_marked_cells() "
                       "only works with parallel distributed triangulations."));
            }
#  endif
        }



      } // namespace

#endif // DEAL_II_WITH_P4EST



      template <class DoFHandlerType>
      ParallelDistributed<DoFHandlerType>::ParallelDistributed(
        DoFHandlerType &dof_handler)
        : dof_handler(&dof_handler)
      {}



      template <class DoFHandlerType>
      NumberCache
      ParallelDistributed<DoFHandlerType>::distribute_dofs() const
      {
#ifndef DEAL_II_WITH_P4EST
        Assert(false, ExcNotImplemented());
        return NumberCache();
#else
        const unsigned int dim      = DoFHandlerType::dimension;
        const unsigned int spacedim = DoFHandlerType::space_dimension;

        parallel::distributed::Triangulation<dim, spacedim> *triangulation =
          (dynamic_cast<parallel::distributed::Triangulation<dim, spacedim> *>(
            const_cast<::Triangulation<dim, spacedim> *>(
              &dof_handler->get_triangulation())));
        Assert(triangulation != nullptr, ExcInternalError());

        const unsigned int n_cpus =
          Utilities::MPI::n_mpi_processes(triangulation->get_communicator());


        /*
           The following algorithm has a number of stages that are all
           documented in the paper that describes the parallel::distributed
           functionality:

           1/ locally enumerate dofs on locally owned cells
           2/ un-numerate those that are on interfaces with ghost
              cells and that we don't own based on the tie-breaking
              criterion; re-enumerate the remaining ones. the
              end result are that we only enumerate locally owned
              DoFs
           3/ shift indices so that each processor has a unique
              range of indices
           4/ for all locally owned cells that are ghost
              cells somewhere else, send our own DoF indices
              to the appropriate set of other processors
         */

        // --------- Phase 1: enumerate dofs on locally owned cells
        const ::types::global_dof_index n_initial_local_dofs =
          Implementation::distribute_dofs(
            triangulation->locally_owned_subdomain(), *dof_handler);

        // --------- Phase 2: un-numerate dofs on interfaces to ghost cells
        //                    that we don't own; re-enumerate the remaining ones

        // start with the identity permutation of indices
        std::vector<::types::global_dof_index> renumbering(
          n_initial_local_dofs);
        for (::types::global_dof_index i = 0; i < renumbering.size(); ++i)
          renumbering[i] = i;

        {
          std::vector<::types::global_dof_index> local_dof_indices;

          for (auto cell : dof_handler->active_cell_iterators())
            if (cell->is_ghost() && (cell->subdomain_id() <
                                     triangulation->locally_owned_subdomain()))
              {
                // we found a neighboring ghost cell whose subdomain
                // is "stronger" than our own subdomain

                // delete all dofs that live there and that we have
                // previously assigned a number to (i.e. the ones on
                // the interface)
                local_dof_indices.resize(cell->get_fe().dofs_per_cell);
                cell->get_dof_indices(local_dof_indices);
                for (auto &local_dof_index : local_dof_indices)
                  if (local_dof_index != numbers::invalid_dof_index)
                    renumbering[local_dof_index] = numbers::invalid_dof_index;
              }
        }


        // make the remaining indices consecutive
        ::types::global_dof_index n_locally_owned_dofs = 0;
        for (auto &new_index : renumbering)
          if (new_index != numbers::invalid_dof_index)
            new_index = n_locally_owned_dofs++;

        // --------- Phase 3: shift indices so that each processor has a unique
        //                    range of indices
        std::vector<::types::global_dof_index>
          n_locally_owned_dofs_per_processor(n_cpus);

        const int ierr =
          MPI_Allgather(&n_locally_owned_dofs,
                        1,
                        DEAL_II_DOF_INDEX_MPI_TYPE,
                        n_locally_owned_dofs_per_processor.data(),
                        1,
                        DEAL_II_DOF_INDEX_MPI_TYPE,
                        triangulation->get_communicator());
        AssertThrowMPI(ierr);

        const ::types::global_dof_index my_shift =
          std::accumulate(n_locally_owned_dofs_per_processor.begin(),
                          n_locally_owned_dofs_per_processor.begin() +
                            triangulation->locally_owned_subdomain(),
                          static_cast<::types::global_dof_index>(0));
        for (auto &new_index : renumbering)
          if (new_index != numbers::invalid_dof_index)
            new_index += my_shift;

        // now re-enumerate all dofs to this shifted and condensed
        // numbering form.  we renumber some dofs as invalid, so
        // choose the nocheck-version.
        Implementation::renumber_dofs(renumbering,
                                      IndexSet(0),
                                      *dof_handler,
                                      false);

        // now a little bit of housekeeping
        const ::types::global_dof_index n_global_dofs =
          std::accumulate(n_locally_owned_dofs_per_processor.begin(),
                          n_locally_owned_dofs_per_processor.end(),
                          ::types::global_dof_index(0));

        std::vector<IndexSet> locally_owned_dofs_per_processor(
          n_cpus, IndexSet(n_global_dofs));
        {
          ::types::global_dof_index current_shift = 0;
          for (unsigned int i = 0; i < n_cpus; ++i)
            {
              locally_owned_dofs_per_processor[i].add_range(
                current_shift,
                current_shift + n_locally_owned_dofs_per_processor[i]);
              current_shift += n_locally_owned_dofs_per_processor[i];
            }
        }
        NumberCache number_cache(locally_owned_dofs_per_processor,
                                 triangulation->locally_owned_subdomain());
        Assert(number_cache
                   .locally_owned_dofs_per_processor
                     [triangulation->locally_owned_subdomain()]
                   .n_elements() == number_cache.n_locally_owned_dofs,
               ExcInternalError());
        Assert(
          !number_cache
              .locally_owned_dofs_per_processor[triangulation
                                                  ->locally_owned_subdomain()]
              .n_elements() ||
            number_cache
                .locally_owned_dofs_per_processor[triangulation
                                                    ->locally_owned_subdomain()]
                .nth_index_in_set(0) == my_shift,
          ExcInternalError());

        // this ends the phase where we enumerate degrees of freedom on
        // each processor. what is missing is communicating DoF indices
        // on ghost cells

        // --------- Phase 4: for all locally owned cells that are ghost
        //                    cells somewhere else, send our own DoF indices
        //                    to the appropriate set of other processors
        {
          std::vector<bool> user_flags;
          triangulation->save_user_flags(user_flags);
          triangulation->clear_user_flags();

          // figure out which cells are ghost cells on which we have
          // to exchange DoF indices
          const std::map<unsigned int, std::set<::types::subdomain_id>>
            vertices_with_ghost_neighbors =
              triangulation->compute_vertices_with_ghost_neighbors();

          // mark all cells that either have to send data (locally
          // owned cells that are adjacent to ghost neighbors in some
          // way) or receive data (all ghost cells) via the user flags
          for (auto cell : dof_handler->active_cell_iterators())
            if (cell->is_locally_owned())
              {
                for (unsigned int v = 0;
                     v < GeometryInfo<dim>::vertices_per_cell;
                     ++v)
                  if (vertices_with_ghost_neighbors.find(cell->vertex_index(
                        v)) != vertices_with_ghost_neighbors.end())
                    {
                      cell->set_user_flag();
                      break;
                    }
              }
            else if (cell->is_ghost())
              cell->set_user_flag();



          // Send and receive cells. After this, only the local cells
          // are marked, that received new data. This has to be
          // communicated in a second communication step.
          //
          // as explained in the 'distributed' paper, this has to be
          // done twice
          communicate_dof_indices_on_marked_cells(
            *dof_handler,
            vertices_with_ghost_neighbors,
            triangulation->coarse_cell_to_p4est_tree_permutation,
            triangulation->p4est_tree_to_coarse_cell_permutation);

          communicate_dof_indices_on_marked_cells(
            *dof_handler,
            vertices_with_ghost_neighbors,
            triangulation->coarse_cell_to_p4est_tree_permutation,
            triangulation->p4est_tree_to_coarse_cell_permutation);

          // at this point, we must have taken care of the data transfer
          // on all cells we had previously marked. verify this
#  ifdef DEBUG
          for (auto cell : dof_handler->active_cell_iterators())
            Assert(cell->user_flag_set() == false, ExcInternalError());
#  endif

          triangulation->load_user_flags(user_flags);
        }

#  ifdef DEBUG
        // check that we are really done
        {
          std::vector<::types::global_dof_index> local_dof_indices;

          for (auto cell : dof_handler->active_cell_iterators())
            if (!cell->is_artificial())
              {
                local_dof_indices.resize(cell->get_fe().dofs_per_cell);
                cell->get_dof_indices(local_dof_indices);
                if (local_dof_indices.end() !=
                    std::find(local_dof_indices.begin(),
                              local_dof_indices.end(),
                              numbers::invalid_dof_index))
                  {
                    if (cell->is_ghost())
                      {
                        Assert(false,
                               ExcMessage(
                                 "A ghost cell ended up with incomplete "
                                 "DoF index information. This should not "
                                 "have happened!"));
                      }
                    else
                      {
                        Assert(
                          false,
                          ExcMessage(
                            "A locally owned cell ended up with incomplete "
                            "DoF index information. This should not "
                            "have happened!"));
                      }
                  }
              }
        }
#  endif // DEBUG
        return number_cache;
#endif   // DEAL_II_WITH_P4EST
      }



      template <class DoFHandlerType>
      std::vector<NumberCache>
      ParallelDistributed<DoFHandlerType>::distribute_mg_dofs() const
      {
#ifndef DEAL_II_WITH_P4EST
        Assert(false, ExcNotImplemented());
        return std::vector<NumberCache>();
#else
        const unsigned int dim      = DoFHandlerType::dimension;
        const unsigned int spacedim = DoFHandlerType::space_dimension;

        parallel::distributed::Triangulation<dim, spacedim> *triangulation =
          (dynamic_cast<parallel::distributed::Triangulation<dim, spacedim> *>(
            const_cast<::Triangulation<dim, spacedim> *>(
              &dof_handler->get_triangulation())));
        Assert(triangulation != nullptr, ExcInternalError());

        AssertThrow((triangulation->settings &
                     parallel::distributed::Triangulation<dim, spacedim>::
                       construct_multigrid_hierarchy),
                    ExcMessage(
                      "Multigrid DoFs can only be distributed on a parallel "
                      "Triangulation if the flag construct_multigrid_hierarchy "
                      "is set in the constructor."));


        const unsigned int n_cpus =
          Utilities::MPI::n_mpi_processes(triangulation->get_communicator());

        // loop over all levels that exist globally (across all
        // processors), even if the current processor does not in fact
        // have any cells on that level or if the local part of the
        // Triangulation has fewer levels. we need to do this because
        // we need to communicate across all processors on all levels
        const unsigned int       n_levels = triangulation->n_global_levels();
        std::vector<NumberCache> number_caches;
        number_caches.reserve(n_levels);
        for (unsigned int level = 0; level < n_levels; ++level)
          {
            NumberCache level_number_cache;

            //* 1. distribute on own subdomain
            const unsigned int n_initial_local_dofs =
              Implementation::distribute_dofs_on_level(
                triangulation->locally_owned_subdomain(), *dof_handler, level);

            //* 2. iterate over ghostcells and kill dofs that are not
            // owned by us
            std::vector<::types::global_dof_index> renumbering(
              n_initial_local_dofs);
            for (::types::global_dof_index i = 0; i < renumbering.size();
                 ++i)
              renumbering[i] = i;

            if (level < triangulation->n_levels())
              {
                std::vector<::types::global_dof_index> local_dof_indices;

                typename DoFHandlerType::level_cell_iterator
                  cell = dof_handler->begin(level),
                  endc = dof_handler->end(level);

                for (; cell != endc; ++cell)
                  if (cell->level_subdomain_id() !=
                        numbers::artificial_subdomain_id &&
                      (cell->level_subdomain_id() <
                       triangulation->locally_owned_subdomain()))
                    {
                      // we found a neighboring ghost cell whose
                      // subdomain is "stronger" than our own
                      // subdomain

                      // delete all dofs that live there and that we
                      // have previously assigned a number to
                      // (i.e. the ones on the interface)
                      local_dof_indices.resize(cell->get_fe().dofs_per_cell);
                      cell->get_mg_dof_indices(local_dof_indices);
                      for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell;
                           ++i)
                        if (local_dof_indices[i] != numbers::invalid_dof_index)
                          renumbering[local_dof_indices[i]] =
                            numbers::invalid_dof_index;
                    }
              }

            // TODO: make this code simpler with the new constructors of
            // NumberCache make indices consecutive
            level_number_cache.n_locally_owned_dofs = 0;
            for (std::vector<::types::global_dof_index>::iterator it =
                   renumbering.begin();
                 it != renumbering.end();
                 ++it)
              if (*it != numbers::invalid_dof_index)
                *it = level_number_cache.n_locally_owned_dofs++;

            //* 3. communicate local dofcount and shift ids to make
            // them unique
            level_number_cache.n_locally_owned_dofs_per_processor.resize(
              n_cpus);

            int ierr = MPI_Allgather(
              &level_number_cache.n_locally_owned_dofs,
              1,
              DEAL_II_DOF_INDEX_MPI_TYPE,
              &level_number_cache.n_locally_owned_dofs_per_processor[0],
              1,
              DEAL_II_DOF_INDEX_MPI_TYPE,
              triangulation->get_communicator());
            AssertThrowMPI(ierr);

            const ::types::global_dof_index shift = std::accumulate(
              level_number_cache.n_locally_owned_dofs_per_processor.begin(),
              level_number_cache.n_locally_owned_dofs_per_processor.begin() +
                triangulation->locally_owned_subdomain(),
              static_cast<::types::global_dof_index>(0));
            for (std::vector<::types::global_dof_index>::iterator it =
                   renumbering.begin();
                 it != renumbering.end();
                 ++it)
              if (*it != numbers::invalid_dof_index)
                (*it) += shift;

            // now re-enumerate all dofs to this shifted and condensed
            // numbering form.  we renumber some dofs as invalid, so
            // choose the nocheck-version of the function
            //
            // of course there is nothing for us to renumber if the
            // level we are currently dealing with doesn't even exist
            // within the current triangulation, so skip renumbering
            // in that case
            if (level < triangulation->n_levels())
              Implementation::renumber_mg_dofs(
                renumbering, IndexSet(0), *dof_handler, level, false);

            // now a little bit of housekeeping
            level_number_cache.n_global_dofs = std::accumulate(
              level_number_cache.n_locally_owned_dofs_per_processor.begin(),
              level_number_cache.n_locally_owned_dofs_per_processor.end(),
              static_cast<::types::global_dof_index>(0));

            level_number_cache.locally_owned_dofs =
              IndexSet(level_number_cache.n_global_dofs);
            level_number_cache.locally_owned_dofs.add_range(
              shift, shift + level_number_cache.n_locally_owned_dofs);
            level_number_cache.locally_owned_dofs.compress();

            // fill global_dof_indexsets
            level_number_cache.locally_owned_dofs_per_processor.resize(n_cpus);
            {
              ::types::global_dof_index current_shift = 0;
              for (unsigned int i = 0; i < n_cpus; ++i)
                {
                  level_number_cache.locally_owned_dofs_per_processor[i] =
                    IndexSet(level_number_cache.n_global_dofs);
                  level_number_cache.locally_owned_dofs_per_processor[i]
                    .add_range(current_shift,
                               current_shift +
                                 level_number_cache
                                   .n_locally_owned_dofs_per_processor[i]);
                  current_shift +=
                    level_number_cache.n_locally_owned_dofs_per_processor[i];
                }
            }
            Assert(level_number_cache
                       .locally_owned_dofs_per_processor
                         [triangulation->locally_owned_subdomain()]
                       .n_elements() == level_number_cache.n_locally_owned_dofs,
                   ExcInternalError());
            Assert(!level_number_cache
                       .locally_owned_dofs_per_processor
                         [triangulation->locally_owned_subdomain()]
                       .n_elements() ||
                     level_number_cache
                         .locally_owned_dofs_per_processor
                           [triangulation->locally_owned_subdomain()]
                         .nth_index_in_set(0) == shift,
                   ExcInternalError());

            number_caches.emplace_back(level_number_cache);
          }


        //* communicate ghost DoFs
        // We mark all ghost cells by setting the user_flag and then request
        // these cells from the corresponding owners. As this information
        // can be incomplete,
        {
          std::vector<bool> user_flags;
          triangulation->save_user_flags(user_flags);
          triangulation->clear_user_flags();

          // mark all ghost cells for transfer
          {
            typename DoFHandlerType::level_cell_iterator cell,
              endc = dof_handler->end();
            for (cell = dof_handler->begin(); cell != endc; ++cell)
              if (cell->level_subdomain_id() !=
                    ::numbers::artificial_subdomain_id &&
                  !cell->is_locally_owned_on_level())
                cell->set_user_flag();
          }

          // Phase 1. Request all marked cells from corresponding owners. If we
          // managed to get every DoF, remove the user_flag, otherwise we
          // will request them again in the step below.
          communicate_mg_ghost_cells(
            *triangulation,
            *dof_handler,
            triangulation->coarse_cell_to_p4est_tree_permutation,
            triangulation->p4est_tree_to_coarse_cell_permutation);

          // have a barrier so that sends from above and below this
          // place are not mixed up.
          //
          // this is necessary because above we just see if there are
          // messages and then receive them, without discriminating
          // where they come from and whether they were sent in phase
          // 1 or 2 in communicate_mg_ghost_cells() on another
          // processor. the need for a global communication step like
          // this barrier could be avoided by receiving messages
          // specifically from those processors from which we expect
          // messages, and by using different tags for phase 1 and 2,
          // but the cost of a barrier is negligible compared to
          // everything else we do here
          const int ierr = MPI_Barrier(triangulation->get_communicator());
          AssertThrowMPI(ierr);

          // Phase 2, only request the cells that were not completed
          // in Phase 1.
          communicate_mg_ghost_cells(
            *triangulation,
            *dof_handler,
            triangulation->coarse_cell_to_p4est_tree_permutation,
            triangulation->p4est_tree_to_coarse_cell_permutation);

#  ifdef DEBUG
          // make sure we have removed all flags:
          {
            typename DoFHandlerType::level_cell_iterator cell,
              endc = dof_handler->end();
            for (cell = dof_handler->begin(); cell != endc; ++cell)
              if (cell->level_subdomain_id() !=
                    ::numbers::artificial_subdomain_id &&
                  !cell->is_locally_owned_on_level())
                Assert(cell->user_flag_set() == false, ExcInternalError());
          }
#  endif

          triangulation->load_user_flags(user_flags);
        }



#  ifdef DEBUG
        // check that we are really done
        {
          std::vector<::types::global_dof_index> local_dof_indices;
          typename DoFHandlerType::level_cell_iterator cell,
            endc = dof_handler->end();

          for (cell = dof_handler->begin(); cell != endc; ++cell)
            if (cell->level_subdomain_id() !=
                ::numbers::artificial_subdomain_id)
              {
                local_dof_indices.resize(cell->get_fe().dofs_per_cell);
                cell->get_mg_dof_indices(local_dof_indices);
                if (local_dof_indices.end() !=
                    std::find(local_dof_indices.begin(),
                              local_dof_indices.end(),
                              numbers::invalid_dof_index))
                  {
                    Assert(false, ExcMessage("not all DoFs got distributed!"));
                  }
              }
        }
#  endif // DEBUG

        return number_caches;

#endif // DEAL_II_WITH_P4EST
      }


      template <class DoFHandlerType>
      NumberCache
      ParallelDistributed<DoFHandlerType>::renumber_dofs(
        const std::vector<::types::global_dof_index> &new_numbers) const
      {
        (void)new_numbers;

        Assert(new_numbers.size() == dof_handler->n_locally_owned_dofs(),
               ExcInternalError());

#ifndef DEAL_II_WITH_P4EST
        Assert(false, ExcNotImplemented());
        return NumberCache();
#else
        const unsigned int dim      = DoFHandlerType::dimension;
        const unsigned int spacedim = DoFHandlerType::space_dimension;

        parallel::distributed::Triangulation<dim, spacedim> *triangulation =
          (dynamic_cast<parallel::distributed::Triangulation<dim, spacedim> *>(
            const_cast<::Triangulation<dim, spacedim> *>(
              &dof_handler->get_triangulation())));
        Assert(triangulation != nullptr, ExcInternalError());


        // First figure out the new set of locally owned DoF indices.
        // If we own no DoFs, we still need to go through this function,
        // but we can skip this calculation.
        //
        // The IndexSet::add_indices() function is substantially more
        // efficient if the set of indices is already sorted because
        // it can then insert ranges instead of individual elements.
        // consequently, pre-sort the array of new indices
        IndexSet my_locally_owned_new_dof_indices(dof_handler->n_dofs());
        if (dof_handler->n_locally_owned_dofs() > 0)
          {
            std::vector<::types::global_dof_index> new_numbers_sorted =
              new_numbers;
            std::sort(new_numbers_sorted.begin(), new_numbers_sorted.end());

            my_locally_owned_new_dof_indices.add_indices(
              new_numbers_sorted.begin(), new_numbers_sorted.end());
            my_locally_owned_new_dof_indices.compress();

            Assert(my_locally_owned_new_dof_indices.n_elements() ==
                     new_numbers.size(),
                   ExcInternalError());
          }

        // delete all knowledge of DoF indices that are not locally
        // owned. we do so by getting DoF indices on cells, checking
        // whether they are locally owned, if not, setting them to
        // an invalid value, and then setting them again on the current
        // cell
        //
        // DoFs we (i) know about, and (ii) don't own locally must be located
        // either on ghost cells, or on the interface between a locally
        // owned cell and a ghost cell. In any case, it is sufficient
        // to kill them only from the ghost side cell, so loop only over
        // ghost cells
        {
          std::vector<::types::global_dof_index> local_dof_indices;

          for (auto cell : dof_handler->active_cell_iterators())
            if (cell->is_ghost())
              {
                local_dof_indices.resize(cell->get_fe().dofs_per_cell);
                cell->get_dof_indices(local_dof_indices);

                for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell; ++i)
                  // delete a DoF index if it has not already been deleted
                  // (e.g., by visiting a neighboring cell, if it is on the
                  // boundary), and if we don't own it
                  if ((local_dof_indices[i] != numbers::invalid_dof_index) &&
                      (!dof_handler->locally_owned_dofs().is_element(
                        local_dof_indices[i])))
                    local_dof_indices[i] = numbers::invalid_dof_index;

                cell->set_dof_indices(local_dof_indices);
              }
        }


        // renumber. Skip when there is nothing to do because we own no DoF.
        if (dof_handler->locally_owned_dofs().n_elements() > 0)
          Implementation::renumber_dofs(new_numbers,
                                        dof_handler->locally_owned_dofs(),
                                        *dof_handler,
                                        false);

        // communicate newly assigned DoF indices to other processors
        // and get the same information for our own ghost cells.
        //
        // this is the same as phase 4 in the distribute_dofs() algorithm
        {
          std::vector<bool> user_flags;
          triangulation->save_user_flags(user_flags);
          triangulation->clear_user_flags();

          // mark all own cells for transfer
          for (auto cell : dof_handler->active_cell_iterators())
            if (!cell->is_artificial())
              cell->set_user_flag();

          // figure out which cells are ghost cells on which we have
          // to exchange DoF indices
          const std::map<unsigned int, std::set<::types::subdomain_id>>
            vertices_with_ghost_neighbors =
              triangulation->compute_vertices_with_ghost_neighbors();


          // Send and receive cells. After this, only the local cells
          // are marked, that received new data. This has to be
          // communicated in a second communication step.
          //
          // as explained in the 'distributed' paper, this has to be
          // done twice
          communicate_dof_indices_on_marked_cells(
            *dof_handler,
            vertices_with_ghost_neighbors,
            triangulation->coarse_cell_to_p4est_tree_permutation,
            triangulation->p4est_tree_to_coarse_cell_permutation);

          communicate_dof_indices_on_marked_cells(
            *dof_handler,
            vertices_with_ghost_neighbors,
            triangulation->coarse_cell_to_p4est_tree_permutation,
            triangulation->p4est_tree_to_coarse_cell_permutation);

          triangulation->load_user_flags(user_flags);
        }

        // the last step is to update the NumberCache, including knowing which
        // processor owns which DoF index. this requires communication.
        //
        // this step is substantially more complicated than it is in
        // distribute_dofs() because the IndexSets of locally owned DoFs
        // after renumbering may not be contiguous any more. for
        // distribute_dofs() it was enough to exchange the starting
        // indices for each processor and the global number of DoFs,
        // but here we actually have to serialize the IndexSet
        // objects and shop them across the network.
        const unsigned int n_cpus =
          Utilities::MPI::n_mpi_processes(triangulation->get_communicator());
        std::vector<IndexSet> locally_owned_dofs_per_processor(
          n_cpus, IndexSet(dof_handler->n_dofs()));
        {
          // serialize our own IndexSet
          std::vector<char> my_data;
          {
#  ifdef DEAL_II_WITH_ZLIB

            boost::iostreams::filtering_ostream out;
            out.push(
              boost::iostreams::gzip_compressor(boost::iostreams::gzip_params(
                boost::iostreams::gzip::best_compression)));
            out.push(boost::iostreams::back_inserter(my_data));

            boost::archive::binary_oarchive archive(out);

            archive << my_locally_owned_new_dof_indices;
            out.flush();
#  else
            std::ostringstream              out;
            boost::archive::binary_oarchive archive(out);
            archive << my_locally_owned_new_dof_indices;
            const std::string &s = out.str();
            my_data.reserve(s.size());
            my_data.assign(s.begin(), s.end());
#  endif
          }

          // determine maximum size of IndexSet
          const unsigned int max_size =
            Utilities::MPI::max(my_data.size(),
                                triangulation->get_communicator());

          // as the MPI_Allgather call will be reading max_size elements, and
          // as this may be past the end of my_data, we need to increase the
          // size of the local buffer. This is filled with zeros.
          my_data.resize(max_size);

          std::vector<char> buffer(max_size * n_cpus);
          const int         ierr = MPI_Allgather(my_data.data(),
                                         max_size,
                                         MPI_BYTE,
                                         buffer.data(),
                                         max_size,
                                         MPI_BYTE,
                                         triangulation->get_communicator());
          AssertThrowMPI(ierr);

          for (unsigned int i = 0; i < n_cpus; ++i)
            if (i == Utilities::MPI::this_mpi_process(
                       triangulation->get_communicator()))
              locally_owned_dofs_per_processor[i] =
                my_locally_owned_new_dof_indices;
            else
              {
                // copy the data previously received into a stringstream
                // object and then read the IndexSet from it
                std::string decompressed_buffer;

                // first decompress the buffer
                {
#  ifdef DEAL_II_WITH_ZLIB

                  boost::iostreams::filtering_ostream decompressing_stream;
                  decompressing_stream.push(
                    boost::iostreams::gzip_decompressor());
                  decompressing_stream.push(
                    boost::iostreams::back_inserter(decompressed_buffer));

                  decompressing_stream.write(&buffer[i * max_size], max_size);
#  else
                  decompressed_buffer.assign(&buffer[i * max_size], max_size);
#  endif
                }

                // then restore the object from the buffer
                std::istringstream              in(decompressed_buffer);
                boost::archive::binary_iarchive archive(in);

                archive >> locally_owned_dofs_per_processor[i];
              }
        }

        return NumberCache(locally_owned_dofs_per_processor,
                           Utilities::MPI::this_mpi_process(
                             triangulation->get_communicator()));
#endif
      }



      template <class DoFHandlerType>
      NumberCache
      ParallelDistributed<DoFHandlerType>::renumber_mg_dofs(
        const unsigned int                          level,
        const std::vector<types::global_dof_index> &new_numbers) const
      {
        // we only implement the case where the multigrid numbers are
        // renumbered within the processor's partition, rather than the most
        // general case
        const std::vector<IndexSet> &index_sets =
          dof_handler->locally_owned_mg_dofs_per_processor(level);

        constexpr int dim      = DoFHandlerType::dimension;
        constexpr int spacedim = DoFHandlerType::space_dimension;
        const parallel::Triangulation<dim, spacedim> *tr =
          (dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
            &this->dof_handler->get_triangulation()));
        Assert(tr != nullptr, ExcInternalError());

#ifdef DEAL_II_WITH_MPI
        const unsigned int my_rank =
          Utilities::MPI::this_mpi_process(tr->get_communicator());

#  ifdef DEBUG
        for (types::global_dof_index i : new_numbers)
          {
            Assert(index_sets[my_rank].is_element(i),
                   ExcNotImplemented(
                     "Renumberings that change the locally owned mg dofs "
                     "partitioning are currently not implemented for "
                     "the multigrid levels"));
          }
#  endif

        // we need to access all locally relevant degrees of freedom. we
        // use Utilities::MPI::Partitioner for handling the data exchange
        // of the new numbers, which is simply the extraction of ghost data
        IndexSet relevant_dofs;
        DoFTools::extract_locally_relevant_level_dofs(*dof_handler,
                                                      level,
                                                      relevant_dofs);
        std::vector<types::global_dof_index> ghosted_new_numbers(
          relevant_dofs.n_elements());
        {
          Utilities::MPI::Partitioner          partitioner(index_sets[my_rank],
                                                  relevant_dofs,
                                                  tr->get_communicator());
          std::vector<types::global_dof_index> temp_array(
            partitioner.n_import_indices());
          const unsigned int       communication_channel = 17;
          std::vector<MPI_Request> requests;
          partitioner.export_to_ghosted_array_start(
            communication_channel,
            make_array_view(new_numbers),
            make_array_view(temp_array),
            ArrayView<types::global_dof_index>(ghosted_new_numbers.data() +
                                                 new_numbers.size(),
                                               partitioner.n_ghost_indices()),
            requests);
          partitioner.export_to_ghosted_array_finish(
            ArrayView<types::global_dof_index>(ghosted_new_numbers.data() +
                                                 new_numbers.size(),
                                               partitioner.n_ghost_indices()),
            requests);

          // we need to fill the indices of the locally owned part into the
          // new numbers array. their right position is somewhere in the
          // middle of the array, so we first copy the ghosted part from
          // smaller ranks to the front, then insert the data in the middle.
          unsigned int n_ghosts_on_smaller_ranks = 0;
          for (std::pair<unsigned int, unsigned int> t :
               partitioner.ghost_targets())
            {
              if (t.first > my_rank)
                break;
              n_ghosts_on_smaller_ranks += t.second;
            }
          if (n_ghosts_on_smaller_ranks > 0)
            {
              Assert(ghosted_new_numbers.data() != nullptr, ExcInternalError());
              std::memmove(ghosted_new_numbers.data(),
                           ghosted_new_numbers.data() + new_numbers.size(),
                           sizeof(types::global_dof_index) *
                             n_ghosts_on_smaller_ranks);
            }
          if (new_numbers.size() > 0)
            {
              Assert(new_numbers.data() != nullptr, ExcInternalError());
              std::memcpy(ghosted_new_numbers.data() +
                            n_ghosts_on_smaller_ranks,
                          new_numbers.data(),
                          sizeof(types::global_dof_index) * new_numbers.size());
            }
        }

        // in case we do not own any of the given level (but only some remote
        // processor), we do not need to call the renumbering
        if (level < this->dof_handler->get_triangulation().n_levels())
          Implementation::renumber_mg_dofs(
            ghosted_new_numbers, relevant_dofs, *dof_handler, level, true);
#else
        (void)new_numbers;
        Assert(false, ExcNotImplemented());
#endif

        return NumberCache(
          index_sets, Utilities::MPI::this_mpi_process(tr->get_communicator()));
      }
    } // namespace Policy
  }   // namespace DoFHandlerImplementation
} // namespace internal



/*-------------- Explicit Instantiations -------------------------------*/
#include "dof_handler_policy.inst"


DEAL_II_NAMESPACE_CLOSE