mg_tools.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 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/logstream.h>
#include <deal.II/base/mg_level_object.h>
#include <deal.II/base/thread_management.h>

#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>

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

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

#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparsity_pattern.h>

#include <deal.II/multigrid/mg_base.h>
#include <deal.II/multigrid/mg_constrained_dofs.h>
#include <deal.II/multigrid/mg_tools.h>

#include <deal.II/numerics/matrix_tools.h>

#include <algorithm>
#include <numeric>
#include <vector>

DEAL_II_NAMESPACE_OPEN


namespace MGTools
{
  // specializations for 1D
  template <>
  void
  compute_row_length_vector(const DoFHandler<1, 1> &,
                            const unsigned int,
                            std::vector<unsigned int> &,
                            const DoFTools::Coupling)
  {
    Assert(false, ExcNotImplemented());
  }



  template <>
  void
  compute_row_length_vector(const DoFHandler<1, 1> &,
                            const unsigned int,
                            std::vector<unsigned int> &,
                            const Table<2, DoFTools::Coupling> &,
                            const Table<2, DoFTools::Coupling> &)
  {
    Assert(false, ExcNotImplemented());
  }



  template <>
  void
  compute_row_length_vector(const DoFHandler<1, 2> &,
                            const unsigned int,
                            std::vector<unsigned int> &,
                            const DoFTools::Coupling)
  {
    Assert(false, ExcNotImplemented());
  }


  template <>
  void
  compute_row_length_vector(const DoFHandler<1, 2> &,
                            const unsigned int,
                            std::vector<unsigned int> &,
                            const Table<2, DoFTools::Coupling> &,
                            const Table<2, DoFTools::Coupling> &)
  {
    Assert(false, ExcNotImplemented());
  }



  // Template for 2D and 3D. For 1D see specialization above
  template <int dim, int spacedim>
  void
  compute_row_length_vector(const DoFHandler<dim, spacedim> &dofs,
                            const unsigned int               level,
                            std::vector<unsigned int> &      row_lengths,
                            const DoFTools::Coupling         flux_coupling)
  {
    Assert(row_lengths.size() == dofs.n_dofs(),
           ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs()));

    // Function starts here by
    // resetting the counters.
    std::fill(row_lengths.begin(), row_lengths.end(), 0);
    // We need the user flags, so we
    // save them for later restoration
    std::vector<bool> old_flags;
    // We need a non-constant
    // triangulation for the user
    // flags. Since we restore them in
    // the end, this cast is safe.
    Triangulation<dim, spacedim> &user_flags_triangulation =
      const_cast<Triangulation<dim, spacedim> &>(dofs.get_triangulation());
    user_flags_triangulation.save_user_flags(old_flags);
    user_flags_triangulation.clear_user_flags();

    const typename DoFHandler<dim, spacedim>::cell_iterator end =
      dofs.end(level);
    typename DoFHandler<dim, spacedim>::active_cell_iterator cell;
    std::vector<types::global_dof_index>                     cell_indices;
    std::vector<types::global_dof_index>                     neighbor_indices;

    // We loop over cells and go from
    // cells to lower dimensional
    // objects. This is the only way to
    // cope with the fact, that an
    // unknown number of cells may
    // share an object of dimension
    // smaller than dim-1.
    for (cell = dofs.begin(level); cell != end; ++cell)
      {
        const FiniteElement<dim> &fe = cell->get_fe();
        cell_indices.resize(fe.dofs_per_cell);
        cell->get_mg_dof_indices(cell_indices);
        unsigned int i = 0;
        // First, dofs on
        // vertices. We assume that
        // each vertex dof couples
        // with all dofs on
        // adjacent grid cells.

        // Adding all dofs of the cells
        // will add dofs of the faces
        // of the cell adjacent to the
        // vertex twice. Therefore, we
        // subtract these here and add
        // them in a loop over the
        // faces below.

        // in 1d, faces and vertices
        // are identical. Nevertheless,
        // this will only work if
        // dofs_per_face is zero and
        // dofs_per_vertex is
        // arbitrary, not the other way
        // round.
        // TODO: This assumes that even in hp context, the dofs per face
        // coincide!
        unsigned int increment = fe.dofs_per_cell - dim * fe.dofs_per_face;
        while (i < fe.first_line_index)
          row_lengths[cell_indices[i++]] += increment;
        // From now on, if an object is
        // a cell, its dofs only couple
        // inside the cell. Since the
        // faces are handled below, we
        // have to subtract ALL faces
        // in this case.

        // In all other cases we
        // subtract adjacent faces to be
        // added in the loop below.
        increment = (dim > 1) ?
                      fe.dofs_per_cell - (dim - 1) * fe.dofs_per_face :
                      fe.dofs_per_cell -
                        GeometryInfo<dim>::faces_per_cell * fe.dofs_per_face;
        while (i < fe.first_quad_index)
          row_lengths[cell_indices[i++]] += increment;

        // Now quads in 2D and 3D
        increment = (dim > 2) ?
                      fe.dofs_per_cell - (dim - 2) * fe.dofs_per_face :
                      fe.dofs_per_cell -
                        GeometryInfo<dim>::faces_per_cell * fe.dofs_per_face;
        while (i < fe.first_hex_index)
          row_lengths[cell_indices[i++]] += increment;
        // Finally, cells in 3D
        increment = fe.dofs_per_cell -
                    GeometryInfo<dim>::faces_per_cell * fe.dofs_per_face;
        while (i < fe.dofs_per_cell)
          row_lengths[cell_indices[i++]] += increment;

        // At this point, we have
        // counted all dofs
        // contributiong from cells
        // coupled topologically to the
        // adjacent cells, but we
        // subtracted some faces.

        // Now, let's go by the faces
        // and add the missing
        // contribution as well as the
        // flux contributions.
        for (unsigned int iface = 0; iface < GeometryInfo<dim>::faces_per_cell;
             ++iface)
          {
            bool level_boundary = cell->at_boundary(iface);
            typename DoFHandler<dim, spacedim>::cell_iterator neighbor;
            if (!level_boundary)
              {
                neighbor = cell->neighbor(iface);
                if (static_cast<unsigned int>(neighbor->level()) != level)
                  level_boundary = true;
              }

            if (level_boundary)
              {
                for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
                     ++local_dof)
                  row_lengths[cell_indices[local_dof]] += fe.dofs_per_face;
                continue;
              }

            const FiniteElement<dim> &nfe = neighbor->get_fe();
            typename DoFHandler<dim, spacedim>::face_iterator face =
              cell->face(iface);

            // Flux couplings are
            // computed from both sides
            // for simplicity.

            // The dofs on the common face
            // will be handled below,
            // therefore, we subtract them
            // here.
            if (flux_coupling != DoFTools::none)
              {
                const unsigned int dof_increment =
                  nfe.dofs_per_cell - nfe.dofs_per_face;
                for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
                     ++local_dof)
                  row_lengths[cell_indices[local_dof]] += dof_increment;
              }

            // Do this only once per
            // face.
            if (face->user_flag_set())
              continue;
            face->set_user_flag();
            // At this point, we assume
            // that each cell added its
            // dofs minus the face to
            // the couplings of the
            // face dofs. Since we
            // subtracted two faces, we
            // have to re-add one.

            // If one side of the face
            // is refined, all the fine
            // face dofs couple with
            // the coarse one.
            neighbor_indices.resize(nfe.dofs_per_cell);
            neighbor->get_mg_dof_indices(neighbor_indices);
            for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
                 ++local_dof)
              row_lengths[cell_indices[local_dof]] += nfe.dofs_per_face;
            for (unsigned int local_dof = 0; local_dof < nfe.dofs_per_cell;
                 ++local_dof)
              row_lengths[neighbor_indices[local_dof]] += fe.dofs_per_face;
          }
      }
    user_flags_triangulation.load_user_flags(old_flags);
  }


  // This is the template for 2D and 3D. See version for 1D above
  template <int dim, int spacedim>
  void
  compute_row_length_vector(const DoFHandler<dim, spacedim> &   dofs,
                            const unsigned int                  level,
                            std::vector<unsigned int> &         row_lengths,
                            const Table<2, DoFTools::Coupling> &couplings,
                            const Table<2, DoFTools::Coupling> &flux_couplings)
  {
    Assert(row_lengths.size() == dofs.n_dofs(),
           ExcDimensionMismatch(row_lengths.size(), dofs.n_dofs()));

    // Function starts here by
    // resetting the counters.
    std::fill(row_lengths.begin(), row_lengths.end(), 0);
    // We need the user flags, so we
    // save them for later restoration
    std::vector<bool> old_flags;
    // We need a non-constant
    // triangulation for the user
    // flags. Since we restore them in
    // the end, this cast is safe.
    Triangulation<dim, spacedim> &user_flags_triangulation =
      const_cast<Triangulation<dim, spacedim> &>(dofs.get_triangulation());
    user_flags_triangulation.save_user_flags(old_flags);
    user_flags_triangulation.clear_user_flags();

    const typename DoFHandler<dim, spacedim>::cell_iterator end =
      dofs.end(level);
    typename DoFHandler<dim, spacedim>::active_cell_iterator cell;
    std::vector<types::global_dof_index>                     cell_indices;
    std::vector<types::global_dof_index>                     neighbor_indices;

    // We have to translate the
    // couplings from components to
    // blocks, so this works for
    // nonprimitive elements as well.
    std::vector<Table<2, DoFTools::Coupling>> couple_cell;
    std::vector<Table<2, DoFTools::Coupling>> couple_face;
    DoFTools::convert_couplings_to_blocks(dofs, couplings, couple_cell);
    DoFTools::convert_couplings_to_blocks(dofs, flux_couplings, couple_face);

    // We loop over cells and go from
    // cells to lower dimensional
    // objects. This is the only way to
    // cope with the fact, that an
    // unknown number of cells may
    // share an object of dimension
    // smaller than dim-1.
    for (cell = dofs.begin_active(); cell != end; ++cell)
      {
        const FiniteElement<dim> &fe       = cell->get_fe();
        const unsigned int        fe_index = cell->active_fe_index();

        Assert(couplings.n_rows() == fe.n_components(),
               ExcDimensionMismatch(couplings.n_rows(), fe.n_components()));
        Assert(couplings.n_cols() == fe.n_components(),
               ExcDimensionMismatch(couplings.n_cols(), fe.n_components()));
        Assert(flux_couplings.n_rows() == fe.n_components(),
               ExcDimensionMismatch(flux_couplings.n_rows(),
                                    fe.n_components()));
        Assert(flux_couplings.n_cols() == fe.n_components(),
               ExcDimensionMismatch(flux_couplings.n_cols(),
                                    fe.n_components()));

        cell_indices.resize(fe.dofs_per_cell);
        cell->get_mg_dof_indices(cell_indices);
        unsigned int i = 0;
        // First, dofs on
        // vertices. We assume that
        // each vertex dof couples
        // with all dofs on
        // adjacent grid cells.

        // Adding all dofs of the cells
        // will add dofs of the faces
        // of the cell adjacent to the
        // vertex twice. Therefore, we
        // subtract these here and add
        // them in a loop over the
        // faces below.

        // in 1d, faces and vertices
        // are identical. Nevertheless,
        // this will only work if
        // dofs_per_face is zero and
        // dofs_per_vertex is
        // arbitrary, not the other way
        // round.
        unsigned int increment;
        while (i < fe.first_line_index)
          {
            for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
              for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
                   ++mult)
                if (couple_cell[fe_index](fe.system_to_block_index(i).first,
                                          fe.first_block_of_base(base) +
                                            mult) != DoFTools::none)
                  {
                    increment = fe.base_element(base).dofs_per_cell -
                                dim * fe.base_element(base).dofs_per_face;
                    row_lengths[cell_indices[i]] += increment;
                  }
            ++i;
          }
        // From now on, if an object is
        // a cell, its dofs only couple
        // inside the cell. Since the
        // faces are handled below, we
        // have to subtract ALL faces
        // in this case.

        // In all other cases we
        // subtract adjacent faces to be
        // added in the loop below.
        while (i < fe.first_quad_index)
          {
            for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
              for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
                   ++mult)
                if (couple_cell[fe_index](fe.system_to_block_index(i).first,
                                          fe.first_block_of_base(base) +
                                            mult) != DoFTools::none)
                  {
                    increment =
                      fe.base_element(base).dofs_per_cell -
                      ((dim > 1) ? (dim - 1) :
                                   GeometryInfo<dim>::faces_per_cell) *
                        fe.base_element(base).dofs_per_face;
                    row_lengths[cell_indices[i]] += increment;
                  }
            ++i;
          }

        // Now quads in 2D and 3D
        while (i < fe.first_hex_index)
          {
            for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
              for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
                   ++mult)
                if (couple_cell[fe_index](fe.system_to_block_index(i).first,
                                          fe.first_block_of_base(base) +
                                            mult) != DoFTools::none)
                  {
                    increment =
                      fe.base_element(base).dofs_per_cell -
                      ((dim > 2) ? (dim - 2) :
                                   GeometryInfo<dim>::faces_per_cell) *
                        fe.base_element(base).dofs_per_face;
                    row_lengths[cell_indices[i]] += increment;
                  }
            ++i;
          }

        // Finally, cells in 3D
        while (i < fe.dofs_per_cell)
          {
            for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
              for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
                   ++mult)
                if (couple_cell[fe_index](fe.system_to_block_index(i).first,
                                          fe.first_block_of_base(base) +
                                            mult) != DoFTools::none)
                  {
                    increment = fe.base_element(base).dofs_per_cell -
                                GeometryInfo<dim>::faces_per_cell *
                                  fe.base_element(base).dofs_per_face;
                    row_lengths[cell_indices[i]] += increment;
                  }
            ++i;
          }

        // At this point, we have
        // counted all dofs
        // contributiong from cells
        // coupled topologically to the
        // adjacent cells, but we
        // subtracted some faces.

        // Now, let's go by the faces
        // and add the missing
        // contribution as well as the
        // flux contributions.
        for (unsigned int iface = 0; iface < GeometryInfo<dim>::faces_per_cell;
             ++iface)
          {
            bool level_boundary = cell->at_boundary(iface);
            typename DoFHandler<dim, spacedim>::cell_iterator neighbor;
            if (!level_boundary)
              {
                neighbor = cell->neighbor(iface);
                if (static_cast<unsigned int>(neighbor->level()) != level)
                  level_boundary = true;
              }

            if (level_boundary)
              {
                for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
                     ++local_dof)
                  row_lengths[cell_indices[local_dof]] += fe.dofs_per_face;
                continue;
              }

            const FiniteElement<dim> &nfe = neighbor->get_fe();
            typename DoFHandler<dim, spacedim>::face_iterator face =
              cell->face(iface);

            // Flux couplings are
            // computed from both sides
            // for simplicity.

            // The dofs on the common face
            // will be handled below,
            // therefore, we subtract them
            // here.
            for (unsigned int base = 0; base < nfe.n_base_elements(); ++base)
              for (unsigned int mult = 0; mult < nfe.element_multiplicity(base);
                   ++mult)
                for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
                     ++local_dof)
                  if (couple_face[fe_index](
                        fe.system_to_block_index(local_dof).first,
                        nfe.first_block_of_base(base) + mult) != DoFTools::none)
                    {
                      const unsigned int dof_increment =
                        nfe.base_element(base).dofs_per_cell -
                        nfe.base_element(base).dofs_per_face;
                      row_lengths[cell_indices[local_dof]] += dof_increment;
                    }

            // Do this only once per
            // face and not on the
            // hanging faces.
            if (face->user_flag_set())
              continue;
            face->set_user_flag();
            // At this point, we assume
            // that each cell added its
            // dofs minus the face to
            // the couplings of the
            // face dofs. Since we
            // subtracted two faces, we
            // have to re-add one.

            // If one side of the face
            // is refined, all the fine
            // face dofs couple with
            // the coarse one.

            // Wolfgang, do they couple
            // with each other by
            // constraints?

            // This will not work with
            // different couplings on
            // different cells.
            neighbor_indices.resize(nfe.dofs_per_cell);
            neighbor->get_mg_dof_indices(neighbor_indices);
            for (unsigned int base = 0; base < nfe.n_base_elements(); ++base)
              for (unsigned int mult = 0; mult < nfe.element_multiplicity(base);
                   ++mult)
                for (unsigned int local_dof = 0; local_dof < fe.dofs_per_cell;
                     ++local_dof)
                  if (couple_cell[fe_index](
                        fe.system_to_component_index(local_dof).first,
                        nfe.first_block_of_base(base) + mult) != DoFTools::none)
                    row_lengths[cell_indices[local_dof]] +=
                      nfe.base_element(base).dofs_per_face;
            for (unsigned int base = 0; base < fe.n_base_elements(); ++base)
              for (unsigned int mult = 0; mult < fe.element_multiplicity(base);
                   ++mult)
                for (unsigned int local_dof = 0; local_dof < nfe.dofs_per_cell;
                     ++local_dof)
                  if (couple_cell[fe_index](
                        nfe.system_to_component_index(local_dof).first,
                        fe.first_block_of_base(base) + mult) != DoFTools::none)
                    row_lengths[neighbor_indices[local_dof]] +=
                      fe.base_element(base).dofs_per_face;
          }
      }
    user_flags_triangulation.load_user_flags(old_flags);
  }



  template <typename DoFHandlerType, typename SparsityPatternType>
  void
  make_sparsity_pattern(const DoFHandlerType &dof,
                        SparsityPatternType & sparsity,
                        const unsigned int    level)
  {
    const types::global_dof_index n_dofs = dof.n_dofs(level);
    (void)n_dofs;

    Assert(sparsity.n_rows() == n_dofs,
           ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
    Assert(sparsity.n_cols() == n_dofs,
           ExcDimensionMismatch(sparsity.n_cols(), n_dofs));

    const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
    std::vector<types::global_dof_index>   dofs_on_this_cell(dofs_per_cell);
    typename DoFHandlerType::cell_iterator cell = dof.begin(level),
                                           endc = dof.end(level);
    for (; cell != endc; ++cell)
      if (dof.get_triangulation().locally_owned_subdomain() ==
            numbers::invalid_subdomain_id ||
          cell->level_subdomain_id() ==
            dof.get_triangulation().locally_owned_subdomain())
        {
          cell->get_mg_dof_indices(dofs_on_this_cell);
          // make sparsity pattern for this cell
          for (unsigned int i = 0; i < dofs_per_cell; ++i)
            for (unsigned int j = 0; j < dofs_per_cell; ++j)
              sparsity.add(dofs_on_this_cell[i], dofs_on_this_cell[j]);
        }
  }



  template <int dim, typename SparsityPatternType, int spacedim>
  void
  make_flux_sparsity_pattern(const DoFHandler<dim, spacedim> &dof,
                             SparsityPatternType &            sparsity,
                             const unsigned int               level)
  {
    const types::global_dof_index n_dofs = dof.n_dofs(level);
    (void)n_dofs;

    Assert(sparsity.n_rows() == n_dofs,
           ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
    Assert(sparsity.n_cols() == n_dofs,
           ExcDimensionMismatch(sparsity.n_cols(), n_dofs));

    const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
    std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
    std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
    typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(level),
                                                      endc = dof.end(level);
    for (; cell != endc; ++cell)
      {
        if (!cell->is_locally_owned_on_level())
          continue;

        cell->get_mg_dof_indices(dofs_on_this_cell);
        // make sparsity pattern for this cell
        for (unsigned int i = 0; i < dofs_per_cell; ++i)
          for (unsigned int j = 0; j < dofs_per_cell; ++j)
            sparsity.add(dofs_on_this_cell[i], dofs_on_this_cell[j]);

        // Loop over all interior neighbors
        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
             ++face)
          {
            bool use_face = false;
            if ((!cell->at_boundary(face)) &&
                (static_cast<unsigned int>(cell->neighbor_level(face)) ==
                 level))
              use_face = true;
            else if (cell->has_periodic_neighbor(face) &&
                     (static_cast<unsigned int>(
                        cell->periodic_neighbor_level(face)) == level))
              use_face = true;

            if (use_face)
              {
                typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
                  cell->neighbor_or_periodic_neighbor(face);
                neighbor->get_mg_dof_indices(dofs_on_other_cell);
                // only add one direction The other is taken care of by
                // neighbor (except when the neighbor is not owned by the same
                // processor)
                for (unsigned int i = 0; i < dofs_per_cell; ++i)
                  {
                    for (unsigned int j = 0; j < dofs_per_cell; ++j)
                      {
                        sparsity.add(dofs_on_this_cell[i],
                                     dofs_on_other_cell[j]);
                      }
                  }
                if (neighbor->is_locally_owned_on_level() == false)
                  for (unsigned int i = 0; i < dofs_per_cell; ++i)
                    for (unsigned int j = 0; j < dofs_per_cell; ++j)
                      {
                        sparsity.add(dofs_on_other_cell[i],
                                     dofs_on_other_cell[j]);
                        sparsity.add(dofs_on_other_cell[i],
                                     dofs_on_this_cell[j]);
                      }
              }
          }
      }
  }



  template <int dim, typename SparsityPatternType, int spacedim>
  void
  make_flux_sparsity_pattern_edge(const DoFHandler<dim, spacedim> &dof,
                                  SparsityPatternType &            sparsity,
                                  const unsigned int               level)
  {
    Assert((level >= 1) && (level < dof.get_triangulation().n_global_levels()),
           ExcIndexRange(level, 1, dof.get_triangulation().n_global_levels()));

    const types::global_dof_index fine_dofs   = dof.n_dofs(level);
    const types::global_dof_index coarse_dofs = dof.n_dofs(level - 1);
    (void)fine_dofs;
    (void)coarse_dofs;

    // Matrix maps from fine level to coarse level

    Assert(sparsity.n_rows() == coarse_dofs,
           ExcDimensionMismatch(sparsity.n_rows(), coarse_dofs));
    Assert(sparsity.n_cols() == fine_dofs,
           ExcDimensionMismatch(sparsity.n_cols(), fine_dofs));

    const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
    std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
    std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
    typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(level),
                                                      endc = dof.end(level);
    for (; cell != endc; ++cell)
      {
        if (!cell->is_locally_owned_on_level())
          continue;

        cell->get_mg_dof_indices(dofs_on_this_cell);
        // Loop over all interior neighbors
        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
             ++face)
          {
            // Neighbor is coarser
            bool use_face = false;
            if ((!cell->at_boundary(face)) &&
                (static_cast<unsigned int>(cell->neighbor_level(face)) !=
                 level))
              use_face = true;
            else if (cell->has_periodic_neighbor(face) &&
                     (static_cast<unsigned int>(
                        cell->periodic_neighbor_level(face)) != level))
              use_face = true;

            if (use_face)
              {
                typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
                  cell->neighbor_or_periodic_neighbor(face);
                neighbor->get_mg_dof_indices(dofs_on_other_cell);

                for (unsigned int i = 0; i < dofs_per_cell; ++i)
                  {
                    for (unsigned int j = 0; j < dofs_per_cell; ++j)
                      {
                        sparsity.add(dofs_on_other_cell[i],
                                     dofs_on_this_cell[j]);
                        sparsity.add(dofs_on_other_cell[j],
                                     dofs_on_this_cell[i]);
                      }
                  }
              }
          }
      }
  }



  template <int dim, typename SparsityPatternType, int spacedim>
  void
  make_flux_sparsity_pattern(const DoFHandler<dim, spacedim> &   dof,
                             SparsityPatternType &               sparsity,
                             const unsigned int                  level,
                             const Table<2, DoFTools::Coupling> &int_mask,
                             const Table<2, DoFTools::Coupling> &flux_mask)
  {
    const FiniteElement<dim> &    fe     = dof.get_fe();
    const types::global_dof_index n_dofs = dof.n_dofs(level);
    const unsigned int            n_comp = fe.n_components();
    (void)n_dofs;
    (void)n_comp;

    Assert(sparsity.n_rows() == n_dofs,
           ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
    Assert(sparsity.n_cols() == n_dofs,
           ExcDimensionMismatch(sparsity.n_cols(), n_dofs));
    Assert(int_mask.n_rows() == n_comp,
           ExcDimensionMismatch(int_mask.n_rows(), n_comp));
    Assert(int_mask.n_cols() == n_comp,
           ExcDimensionMismatch(int_mask.n_cols(), n_comp));
    Assert(flux_mask.n_rows() == n_comp,
           ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
    Assert(flux_mask.n_cols() == n_comp,
           ExcDimensionMismatch(flux_mask.n_cols(), n_comp));

    const unsigned int                   total_dofs = fe.dofs_per_cell;
    std::vector<types::global_dof_index> dofs_on_this_cell(total_dofs);
    std::vector<types::global_dof_index> dofs_on_other_cell(total_dofs);
    Table<2, bool>                       support_on_face(total_dofs,
                                   GeometryInfo<dim>::faces_per_cell);

    typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(level),
                                                      endc = dof.end(level);

    const Table<2, DoFTools::Coupling>
      int_dof_mask =
        DoFTools::dof_couplings_from_component_couplings(fe, int_mask),
      flux_dof_mask =
        DoFTools::dof_couplings_from_component_couplings(fe, flux_mask);

    for (unsigned int i = 0; i < total_dofs; ++i)
      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
        support_on_face(i, f) = fe.has_support_on_face(i, f);

    // 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;
    dof.get_triangulation().save_user_flags(user_flags);
    const_cast<Triangulation<dim, spacedim> &>(dof.get_triangulation())
      .clear_user_flags();

    for (; cell != endc; ++cell)
      {
        if (!cell->is_locally_owned_on_level())
          continue;

        cell->get_mg_dof_indices(dofs_on_this_cell);
        // make sparsity pattern for this cell
        for (unsigned int i = 0; i < total_dofs; ++i)
          for (unsigned int j = 0; j < total_dofs; ++j)
            if (int_dof_mask[i][j] != DoFTools::none)
              sparsity.add(dofs_on_this_cell[i], dofs_on_this_cell[j]);

        // Loop over all interior neighbors
        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
             ++face)
          {
            typename DoFHandler<dim, spacedim>::face_iterator cell_face =
              cell->face(face);
            if (cell_face->user_flag_set())
              continue;

            if (cell->at_boundary(face) && !cell->has_periodic_neighbor(face))
              {
                for (unsigned int i = 0; i < total_dofs; ++i)
                  {
                    const bool i_non_zero_i = support_on_face(i, face);
                    for (unsigned int j = 0; j < total_dofs; ++j)
                      {
                        const bool j_non_zero_i = support_on_face(j, face);

                        if (flux_dof_mask(i, j) == DoFTools::always)
                          sparsity.add(dofs_on_this_cell[i],
                                       dofs_on_this_cell[j]);
                        if (flux_dof_mask(i, j) == DoFTools::nonzero &&
                            i_non_zero_i && j_non_zero_i)
                          sparsity.add(dofs_on_this_cell[i],
                                       dofs_on_this_cell[j]);
                      }
                  }
              }
            else
              {
                typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
                  cell->neighbor_or_periodic_neighbor(face);

                if (neighbor->level() < cell->level())
                  continue;

                unsigned int neighbor_face =
                  cell->has_periodic_neighbor(face) ?
                    cell->periodic_neighbor_of_periodic_neighbor(face) :
                    cell->neighbor_of_neighbor(face);

                neighbor->get_mg_dof_indices(dofs_on_other_cell);
                for (unsigned int i = 0; i < total_dofs; ++i)
                  {
                    const bool i_non_zero_i = support_on_face(i, face);
                    const bool i_non_zero_e = support_on_face(i, neighbor_face);
                    for (unsigned int j = 0; j < total_dofs; ++j)
                      {
                        const bool j_non_zero_i = support_on_face(j, face);
                        const bool j_non_zero_e =
                          support_on_face(j, neighbor_face);
                        if (flux_dof_mask(i, j) == DoFTools::always)
                          {
                            sparsity.add(dofs_on_this_cell[i],
                                         dofs_on_other_cell[j]);
                            sparsity.add(dofs_on_other_cell[i],
                                         dofs_on_this_cell[j]);
                            sparsity.add(dofs_on_this_cell[i],
                                         dofs_on_this_cell[j]);
                            sparsity.add(dofs_on_other_cell[i],
                                         dofs_on_other_cell[j]);
                          }
                        if (flux_dof_mask(i, j) == DoFTools::nonzero)
                          {
                            if (i_non_zero_i && j_non_zero_e)
                              sparsity.add(dofs_on_this_cell[i],
                                           dofs_on_other_cell[j]);
                            if (i_non_zero_e && j_non_zero_i)
                              sparsity.add(dofs_on_other_cell[i],
                                           dofs_on_this_cell[j]);
                            if (i_non_zero_i && j_non_zero_i)
                              sparsity.add(dofs_on_this_cell[i],
                                           dofs_on_this_cell[j]);
                            if (i_non_zero_e && j_non_zero_e)
                              sparsity.add(dofs_on_other_cell[i],
                                           dofs_on_other_cell[j]);
                          }

                        if (flux_dof_mask(j, i) == DoFTools::always)
                          {
                            sparsity.add(dofs_on_this_cell[j],
                                         dofs_on_other_cell[i]);
                            sparsity.add(dofs_on_other_cell[j],
                                         dofs_on_this_cell[i]);
                            sparsity.add(dofs_on_this_cell[j],
                                         dofs_on_this_cell[i]);
                            sparsity.add(dofs_on_other_cell[j],
                                         dofs_on_other_cell[i]);
                          }
                        if (flux_dof_mask(j, i) == DoFTools::nonzero)
                          {
                            if (j_non_zero_i && i_non_zero_e)
                              sparsity.add(dofs_on_this_cell[j],
                                           dofs_on_other_cell[i]);
                            if (j_non_zero_e && i_non_zero_i)
                              sparsity.add(dofs_on_other_cell[j],
                                           dofs_on_this_cell[i]);
                            if (j_non_zero_i && i_non_zero_i)
                              sparsity.add(dofs_on_this_cell[j],
                                           dofs_on_this_cell[i]);
                            if (j_non_zero_e && i_non_zero_e)
                              sparsity.add(dofs_on_other_cell[j],
                                           dofs_on_other_cell[i]);
                          }
                      }
                  }
                neighbor->face(neighbor_face)->set_user_flag();
              }
          }
      }

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



  template <int dim, typename SparsityPatternType, int spacedim>
  void
  make_flux_sparsity_pattern_edge(const DoFHandler<dim, spacedim> &   dof,
                                  SparsityPatternType &               sparsity,
                                  const unsigned int                  level,
                                  const Table<2, DoFTools::Coupling> &flux_mask)
  {
    const FiniteElement<dim> &fe     = dof.get_fe();
    const unsigned int        n_comp = fe.n_components();
    (void)n_comp;

    Assert((level >= 1) && (level < dof.get_triangulation().n_global_levels()),
           ExcIndexRange(level, 1, dof.get_triangulation().n_global_levels()));

    const types::global_dof_index fine_dofs   = dof.n_dofs(level);
    const types::global_dof_index coarse_dofs = dof.n_dofs(level - 1);
    (void)fine_dofs;
    (void)coarse_dofs;

    // Matrix maps from fine level to coarse level

    Assert(sparsity.n_rows() == coarse_dofs,
           ExcDimensionMismatch(sparsity.n_rows(), coarse_dofs));
    Assert(sparsity.n_cols() == fine_dofs,
           ExcDimensionMismatch(sparsity.n_cols(), fine_dofs));
    Assert(flux_mask.n_rows() == n_comp,
           ExcDimensionMismatch(flux_mask.n_rows(), n_comp));
    Assert(flux_mask.n_cols() == n_comp,
           ExcDimensionMismatch(flux_mask.n_cols(), n_comp));

    const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
    std::vector<types::global_dof_index> dofs_on_this_cell(dofs_per_cell);
    std::vector<types::global_dof_index> dofs_on_other_cell(dofs_per_cell);
    Table<2, bool>                       support_on_face(dofs_per_cell,
                                   GeometryInfo<dim>::faces_per_cell);

    typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(level),
                                                      endc = dof.end(level);

    const Table<2, DoFTools::Coupling> flux_dof_mask =
      DoFTools::dof_couplings_from_component_couplings(fe, flux_mask);

    for (unsigned int i = 0; i < dofs_per_cell; ++i)
      for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
        support_on_face(i, f) = fe.has_support_on_face(i, f);

    for (; cell != endc; ++cell)
      {
        if (!cell->is_locally_owned_on_level())
          continue;

        cell->get_mg_dof_indices(dofs_on_this_cell);
        // Loop over all interior neighbors
        for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
             ++face)
          {
            // Neighbor is coarser
            bool use_face = false;
            if ((!cell->at_boundary(face)) &&
                (static_cast<unsigned int>(cell->neighbor_level(face)) !=
                 level))
              use_face = true;
            else if (cell->has_periodic_neighbor(face) &&
                     (static_cast<unsigned int>(
                        cell->periodic_neighbor_level(face)) != level))
              use_face = true;

            if (use_face)
              {
                typename DoFHandler<dim, spacedim>::cell_iterator neighbor =
                  cell->neighbor_or_periodic_neighbor(face);
                neighbor->get_mg_dof_indices(dofs_on_other_cell);

                for (unsigned int i = 0; i < dofs_per_cell; ++i)
                  {
                    for (unsigned int j = 0; j < dofs_per_cell; ++j)
                      {
                        if (flux_dof_mask(i, j) != DoFTools::none)
                          {
                            sparsity.add(dofs_on_other_cell[i],
                                         dofs_on_this_cell[j]);
                            sparsity.add(dofs_on_other_cell[j],
                                         dofs_on_this_cell[i]);
                          }
                      }
                  }
              }
          }
      }
  }



  template <typename DoFHandlerType, typename SparsityPatternType>
  void
  make_interface_sparsity_pattern(const DoFHandlerType &   dof,
                                  const MGConstrainedDoFs &mg_constrained_dofs,
                                  SparsityPatternType &    sparsity,
                                  const unsigned int       level)
  {
    const types::global_dof_index n_dofs = dof.n_dofs(level);
    (void)n_dofs;

    Assert(sparsity.n_rows() == n_dofs,
           ExcDimensionMismatch(sparsity.n_rows(), n_dofs));
    Assert(sparsity.n_cols() == n_dofs,
           ExcDimensionMismatch(sparsity.n_cols(), n_dofs));

    const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
    std::vector<types::global_dof_index>   dofs_on_this_cell(dofs_per_cell);
    typename DoFHandlerType::cell_iterator cell = dof.begin(level),
                                           endc = dof.end(level);
    for (; cell != endc; ++cell)
      if (cell->is_locally_owned_on_level())
        {
          cell->get_mg_dof_indices(dofs_on_this_cell);
          for (unsigned int i = 0; i < dofs_per_cell; ++i)
            for (unsigned int j = 0; j < dofs_per_cell; ++j)
              if (mg_constrained_dofs.is_interface_matrix_entry(
                    level, dofs_on_this_cell[i], dofs_on_this_cell[j]))
                sparsity.add(dofs_on_this_cell[i], dofs_on_this_cell[j]);
        }
  }



  template <int dim, int spacedim>
  void
  count_dofs_per_component(
    const DoFHandler<dim, spacedim> &                  dof_handler,
    std::vector<std::vector<types::global_dof_index>> &result,
    bool                                               only_once,
    std::vector<unsigned int>                          target_component)
  {
    const FiniteElement<dim> &fe           = dof_handler.get_fe();
    const unsigned int        n_components = fe.n_components();
    const unsigned int        nlevels =
      dof_handler.get_triangulation().n_global_levels();

    Assert(result.size() == nlevels,
           ExcDimensionMismatch(result.size(), nlevels));

    if (target_component.size() == 0)
      {
        target_component.resize(n_components);
        for (unsigned int i = 0; i < n_components; ++i)
          target_component[i] = i;
      }

    Assert(target_component.size() == n_components,
           ExcDimensionMismatch(target_component.size(), n_components));

    for (unsigned int l = 0; l < nlevels; ++l)
      {
        result[l].resize(n_components);
        std::fill(result[l].begin(), result[l].end(), 0U);

        // special case for only one
        // component. treat this first
        // since it does not require any
        // computations
        if (n_components == 1)
          {
            result[l][0] = dof_handler.n_dofs(l);
          }
        else
          {
            // otherwise determine the number
            // of dofs in each component
            // separately. do so in parallel
            std::vector<std::vector<bool>> dofs_in_component(
              n_components, std::vector<bool>(dof_handler.n_dofs(l), false));
            std::vector<ComponentMask> component_select(n_components);
            Threads::TaskGroup<>       tasks;
            for (unsigned int i = 0; i < n_components; ++i)
              {
                void (*fun_ptr)(const unsigned int level,
                                const DoFHandler<dim, spacedim> &,
                                const ComponentMask &,
                                std::vector<bool> &) =
                  &DoFTools::extract_level_dofs<DoFHandler<dim, spacedim>>;

                std::vector<bool> tmp(n_components, false);
                tmp[i]              = true;
                component_select[i] = ComponentMask(tmp);

                tasks += Threads::new_task(fun_ptr,
                                           l,
                                           dof_handler,
                                           component_select[i],
                                           dofs_in_component[i]);
              }
            tasks.join_all();

            // next count what we got
            unsigned int component = 0;
            for (unsigned int b = 0; b < fe.n_base_elements(); ++b)
              {
                const FiniteElement<dim> &base = fe.base_element(b);
                // Dimension of base element
                unsigned int d = base.n_components();

                for (unsigned int m = 0; m < fe.element_multiplicity(b); ++m)
                  {
                    for (unsigned int dd = 0; dd < d; ++dd)
                      {
                        if (base.is_primitive() || (!only_once || dd == 0))
                          result[l][target_component[component]] +=
                            std::count(dofs_in_component[component].begin(),
                                       dofs_in_component[component].end(),
                                       true);
                        ++component;
                      }
                  }
              }
            // finally sanity check
            Assert(!dof_handler.get_fe().is_primitive() ||
                     std::accumulate(result[l].begin(),
                                     result[l].end(),
                                     types::global_dof_index(0)) ==
                       dof_handler.n_dofs(l),
                   ExcInternalError());
          }
      }
  }



  template <typename DoFHandlerType>
  void
  count_dofs_per_block(
    const DoFHandlerType &                             dof_handler,
    std::vector<std::vector<types::global_dof_index>> &dofs_per_block,
    std::vector<unsigned int>                          target_block)
  {
    const FiniteElement<DoFHandlerType::dimension,
                        DoFHandlerType::space_dimension> &fe =
      dof_handler.get_fe();
    const unsigned int n_blocks = fe.n_blocks();
    const unsigned int n_levels =
      dof_handler.get_triangulation().n_global_levels();

    AssertDimension(dofs_per_block.size(), n_levels);

    for (unsigned int l = 0; l < n_levels; ++l)
      std::fill(dofs_per_block[l].begin(), dofs_per_block[l].end(), 0U);
    // If the empty vector was given as
    // default argument, set up this
    // vector as identity.
    if (target_block.size() == 0)
      {
        target_block.resize(n_blocks);
        for (unsigned int i = 0; i < n_blocks; ++i)
          target_block[i] = i;
      }
    Assert(target_block.size() == n_blocks,
           ExcDimensionMismatch(target_block.size(), n_blocks));

    const unsigned int max_block =
      *std::max_element(target_block.begin(), target_block.end());
    const unsigned int n_target_blocks = max_block + 1;
    (void)n_target_blocks;

    for (unsigned int l = 0; l < n_levels; ++l)
      AssertDimension(dofs_per_block[l].size(), n_target_blocks);

    // special case for only one
    // block. treat this first
    // since it does not require any
    // computations
    if (n_blocks == 1)
      {
        for (unsigned int l = 0; l < n_levels; ++l)
          dofs_per_block[l][0] = dof_handler.n_dofs(l);
        return;
      }
    // otherwise determine the number
    // of dofs in each block
    // separately. do so in parallel
    for (unsigned int l = 0; l < n_levels; ++l)
      {
        std::vector<std::vector<bool>> dofs_in_block(
          n_blocks, std::vector<bool>(dof_handler.n_dofs(l), false));
        std::vector<BlockMask> block_select(n_blocks);
        Threads::TaskGroup<>   tasks;
        for (unsigned int i = 0; i < n_blocks; ++i)
          {
            void (*fun_ptr)(const unsigned int level,
                            const DoFHandlerType &,
                            const BlockMask &,
                            std::vector<bool> &) =
              &DoFTools::extract_level_dofs<DoFHandlerType>;

            std::vector<bool> tmp(n_blocks, false);
            tmp[i]          = true;
            block_select[i] = tmp;

            tasks += Threads::new_task(
              fun_ptr, l, dof_handler, block_select[i], dofs_in_block[i]);
          }
        tasks.join_all();

        // next count what we got
        for (unsigned int block = 0; block < fe.n_blocks(); ++block)
          dofs_per_block[l][target_block[block]] +=
            std::count(dofs_in_block[block].begin(),
                       dofs_in_block[block].end(),
                       true);
      }
  }



  template <int dim, int spacedim>
  void
  make_boundary_list(
    const DoFHandler<dim, spacedim> &dof,
    const std::map<types::boundary_id, const Function<spacedim> *>
      &                                             function_map,
    std::vector<std::set<types::global_dof_index>> &boundary_indices,
    const ComponentMask &                           component_mask)
  {
    Assert(boundary_indices.size() == dof.get_triangulation().n_global_levels(),
           ExcDimensionMismatch(boundary_indices.size(),
                                dof.get_triangulation().n_global_levels()));

    std::set<types::boundary_id> boundary_ids;
    for (const auto &boundary_function : function_map)
      boundary_ids.insert(boundary_function.first);

    std::vector<IndexSet> boundary_indexset;
    make_boundary_list(dof, boundary_ids, boundary_indexset, component_mask);
    for (unsigned int i = 0; i < dof.get_triangulation().n_global_levels(); ++i)
      boundary_indices[i].insert(boundary_indexset[i].begin(),
                                 boundary_indexset[i].end());
  }


  template <int dim, int spacedim>
  void
  make_boundary_list(const DoFHandler<dim, spacedim> &           dof,
                     const std::map<types::boundary_id,
                                    const Function<spacedim> *> &function_map,
                     std::vector<IndexSet> &boundary_indices,
                     const ComponentMask &  component_mask)
  {
    Assert(boundary_indices.size() == dof.get_triangulation().n_global_levels(),
           ExcDimensionMismatch(boundary_indices.size(),
                                dof.get_triangulation().n_global_levels()));

    std::set<types::boundary_id> boundary_ids;
    for (const auto &boundary_function : function_map)
      boundary_ids.insert(boundary_function.first);

    make_boundary_list(dof, boundary_ids, boundary_indices, component_mask);
  }



  template <int dim, int spacedim>
  void
  make_boundary_list(const DoFHandler<dim, spacedim> &   dof,
                     const std::set<types::boundary_id> &boundary_ids,
                     std::vector<IndexSet> &             boundary_indices,
                     const ComponentMask &               component_mask)
  {
    boundary_indices.resize(dof.get_triangulation().n_global_levels());

    // if for whatever reason we were passed an empty set, return immediately
    if (boundary_ids.size() == 0)
      return;

    for (unsigned int i = 0; i < dof.get_triangulation().n_global_levels(); ++i)
      if (boundary_indices[i].size() == 0)
        boundary_indices[i] = IndexSet(dof.n_dofs(i));

    const unsigned int n_components = DoFTools::n_components(dof);
    const bool         fe_is_system = (n_components != 1);

    std::vector<types::global_dof_index> local_dofs;
    local_dofs.reserve(DoFTools::max_dofs_per_face(dof));
    std::fill(local_dofs.begin(), local_dofs.end(), numbers::invalid_dof_index);

    // First, deal with the simpler case when we have to identify all boundary
    // dofs
    if (component_mask.n_selected_components(n_components) == n_components)
      {
        typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(),
                                                          endc = dof.end();
        for (; cell != endc; ++cell)
          {
            if (dof.get_triangulation().locally_owned_subdomain() !=
                  numbers::invalid_subdomain_id &&
                cell->level_subdomain_id() == numbers::artificial_subdomain_id)
              continue;
            const FiniteElement<dim> &fe    = cell->get_fe();
            const unsigned int        level = cell->level();
            local_dofs.resize(fe.dofs_per_face);

            for (unsigned int face_no = 0;
                 face_no < GeometryInfo<dim>::faces_per_cell;
                 ++face_no)
              if (cell->at_boundary(face_no) == true)
                {
                  const typename DoFHandler<dim, spacedim>::face_iterator face =
                    cell->face(face_no);
                  const types::boundary_id bi = face->boundary_id();
                  // Face is listed in boundary map
                  if (boundary_ids.find(bi) != boundary_ids.end())
                    {
                      face->get_mg_dof_indices(level, local_dofs);
                      boundary_indices[level].add_indices(local_dofs.begin(),
                                                          local_dofs.end());
                    }
                }
          }
      }
    else
      {
        Assert(component_mask.n_selected_components(n_components) > 0,
               ExcMessage(
                 "It's probably worthwhile to select at least one component."));

        typename DoFHandler<dim, spacedim>::cell_iterator cell = dof.begin(),
                                                          endc = dof.end();
        for (; cell != endc; ++cell)
          if (dof.get_triangulation().locally_owned_subdomain() ==
                numbers::invalid_subdomain_id ||
              cell->level_subdomain_id() != numbers::artificial_subdomain_id)
            for (unsigned int face_no = 0;
                 face_no < GeometryInfo<dim>::faces_per_cell;
                 ++face_no)
              {
                if (cell->at_boundary(face_no) == false)
                  continue;

                const FiniteElement<dim> &fe    = cell->get_fe();
                const unsigned int        level = cell->level();

                typename DoFHandler<dim, spacedim>::face_iterator face =
                  cell->face(face_no);
                const types::boundary_id boundary_component =
                  face->boundary_id();
                if (boundary_ids.find(boundary_component) != boundary_ids.end())
                  // we want to constrain this boundary
                  {
                    for (unsigned int i = 0; i < cell->get_fe().dofs_per_cell;
                         ++i)
                      {
                        const ComponentMask &nonzero_component_array =
                          cell->get_fe().get_nonzero_components(i);
                        // if we want to constrain one of the nonzero
                        // components, we have to constrain all of them

                        bool selected = false;
                        for (unsigned int c = 0; c < n_components; ++c)
                          if (nonzero_component_array[c] == true &&
                              component_mask[c] == true)
                            {
                              selected = true;
                              break;
                            }
                        if (selected)
                          for (unsigned int c = 0; c < n_components; ++c)
                            Assert(
                              nonzero_component_array[c] == false ||
                                component_mask[c] == true,
                              ExcMessage(
                                "You are using a non-primitive FiniteElement "
                                "and try to constrain just some of its components!"));
                      }

                    // get indices, physical location and boundary values of
                    // dofs on this face
                    local_dofs.resize(fe.dofs_per_face);
                    face->get_mg_dof_indices(level, local_dofs);
                    if (fe_is_system)
                      {
                        for (unsigned int i = 0; i < local_dofs.size(); ++i)
                          {
                            unsigned int component =
                              numbers::invalid_unsigned_int;
                            if (fe.is_primitive())
                              component =
                                fe.face_system_to_component_index(i).first;
                            else
                              {
                                // Just pick the first of the components
                                // We already know that either all or none
                                // of the components are selected
                                const ComponentMask &nonzero_component_array =
                                  cell->get_fe().get_nonzero_components(i);
                                for (unsigned int c = 0; c < n_components; ++c)
                                  if (nonzero_component_array[c] == true)
                                    {
                                      component = c;
                                      break;
                                    }
                              }
                            Assert(component != numbers::invalid_unsigned_int,
                                   ExcInternalError());
                            if (component_mask[component] == true)
                              boundary_indices[level].add_index(local_dofs[i]);
                          }
                      }
                    else
                      boundary_indices[level].add_indices(local_dofs.begin(),
                                                          local_dofs.end());
                  }
              }
      }
  }


  template <int dim, int spacedim>
  void
  extract_non_interface_dofs(
    const DoFHandler<dim, spacedim> &               mg_dof_handler,
    std::vector<std::set<types::global_dof_index>> &non_interface_dofs)
  {
    Assert(non_interface_dofs.size() ==
             mg_dof_handler.get_triangulation().n_global_levels(),
           ExcDimensionMismatch(
             non_interface_dofs.size(),
             mg_dof_handler.get_triangulation().n_global_levels()));

    const FiniteElement<dim, spacedim> &fe = mg_dof_handler.get_fe();

    const unsigned int dofs_per_cell = fe.dofs_per_cell;
    const unsigned int dofs_per_face = fe.dofs_per_face;

    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
    std::vector<bool>                    cell_dofs(dofs_per_cell, false);
    std::vector<bool> cell_dofs_interface(dofs_per_cell, false);

    typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
                                            endc = mg_dof_handler.end();


    for (; cell != endc; ++cell)
      {
        if (mg_dof_handler.get_triangulation().locally_owned_subdomain() !=
              numbers::invalid_subdomain_id &&
            cell->level_subdomain_id() !=
              mg_dof_handler.get_triangulation().locally_owned_subdomain())
          continue;

        std::fill(cell_dofs.begin(), cell_dofs.end(), false);
        std::fill(cell_dofs_interface.begin(),
                  cell_dofs_interface.end(),
                  false);

        for (unsigned int face_nr = 0;
             face_nr < GeometryInfo<dim>::faces_per_cell;
             ++face_nr)
          {
            const typename DoFHandler<dim, spacedim>::face_iterator face =
              cell->face(face_nr);
            if (!face->at_boundary())
              {
                // interior face
                const typename DoFHandler<dim>::cell_iterator neighbor =
                  cell->neighbor(face_nr);

                if ((neighbor->level() < cell->level()))
                  {
                    for (unsigned int j = 0; j < dofs_per_face; ++j)
                      cell_dofs_interface[fe.face_to_cell_index(j, face_nr)] =
                        true;
                  }
                else
                  {
                    for (unsigned int j = 0; j < dofs_per_face; ++j)
                      cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
                  }
              }
            else
              {
                // boundary face
                for (unsigned int j = 0; j < dofs_per_face; ++j)
                  cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;
              }
          }

        const unsigned int level = cell->level();
        cell->get_mg_dof_indices(local_dof_indices);

        for (unsigned int i = 0; i < dofs_per_cell; ++i)
          if (cell_dofs[i] && !cell_dofs_interface[i])
            non_interface_dofs[level].insert(local_dof_indices[i]);
      }
  }


  template <int dim, int spacedim>
  void
  extract_inner_interface_dofs(const DoFHandler<dim, spacedim> &mg_dof_handler,
                               std::vector<IndexSet> &          interface_dofs)
  {
    Assert(interface_dofs.size() ==
             mg_dof_handler.get_triangulation().n_global_levels(),
           ExcDimensionMismatch(
             interface_dofs.size(),
             mg_dof_handler.get_triangulation().n_global_levels()));

    std::vector<std::vector<types::global_dof_index>> tmp_interface_dofs(
      interface_dofs.size());

    const FiniteElement<dim, spacedim> &fe = mg_dof_handler.get_fe();

    const unsigned int dofs_per_cell = fe.dofs_per_cell;
    const unsigned int dofs_per_face = fe.dofs_per_face;

    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);

    std::vector<bool> cell_dofs(dofs_per_cell, false);

    typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
                                            endc = mg_dof_handler.end();

    for (; cell != endc; ++cell)
      {
        // Do not look at artificial level cells (in a serial computation we
        // need to ignore the level_subdomain_id() because it is never set).
        if (mg_dof_handler.get_triangulation().locally_owned_subdomain() !=
              numbers::invalid_subdomain_id &&
            cell->level_subdomain_id() == numbers::artificial_subdomain_id)
          continue;

        bool has_coarser_neighbor = false;

        std::fill(cell_dofs.begin(), cell_dofs.end(), false);

        for (unsigned int face_nr = 0;
             face_nr < GeometryInfo<dim>::faces_per_cell;
             ++face_nr)
          {
            const typename DoFHandler<dim, spacedim>::face_iterator face =
              cell->face(face_nr);
            if (!face->at_boundary())
              {
                // interior face
                const typename DoFHandler<dim>::cell_iterator neighbor =
                  cell->neighbor(face_nr);

                // only process cell pairs if one or both of them are owned by
                // me (ignore if running in serial)
                if (mg_dof_handler.get_triangulation()
                        .locally_owned_subdomain() !=
                      numbers::invalid_subdomain_id &&
                    neighbor->level_subdomain_id() ==
                      numbers::artificial_subdomain_id)
                  continue;

                // Do refinement face from the coarse side
                if (neighbor->level() < cell->level())
                  {
                    for (unsigned int j = 0; j < dofs_per_face; ++j)
                      cell_dofs[fe.face_to_cell_index(j, face_nr)] = true;

                    has_coarser_neighbor = true;
                  }
              }
          }

        if (has_coarser_neighbor == false)
          continue;

        const unsigned int level = cell->level();
        cell->get_mg_dof_indices(local_dof_indices);

        for (unsigned int i = 0; i < dofs_per_cell; ++i)
          {
            if (cell_dofs[i])
              tmp_interface_dofs[level].push_back(local_dof_indices[i]);
          }
      }

    for (unsigned int l = 0;
         l < mg_dof_handler.get_triangulation().n_global_levels();
         ++l)
      {
        interface_dofs[l].clear();
        std::sort(tmp_interface_dofs[l].begin(), tmp_interface_dofs[l].end());
        interface_dofs[l].add_indices(tmp_interface_dofs[l].begin(),
                                      std::unique(tmp_interface_dofs[l].begin(),
                                                  tmp_interface_dofs[l].end()));
        interface_dofs[l].compress();
      }
  }



  template <int dim, int spacedim>
  unsigned int
  max_level_for_coarse_mesh(const Triangulation<dim, spacedim> &tria)
  {
    // Find minimum level for an active cell in
    // this locally owned subdomain
    // Note: with the way active cells are traversed,
    // the first locally owned cell we find will have
    // the lowest level in the particular subdomain.
    unsigned int min_level = tria.n_global_levels();
    typename Triangulation<dim, spacedim>::active_cell_iterator
      cell = tria.begin_active(),
      endc = tria.end();
    for (; cell != endc; ++cell)
      if (cell->is_locally_owned())
        {
          min_level = cell->level();
          break;
        }

    unsigned int global_min = min_level;
    // If necessary, communicate to find minimum
    // level for an active cell over all subdomains
    if (const parallel::Triangulation<dim, spacedim> *tr =
          dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(&tria))
      global_min = Utilities::MPI::min(min_level, tr->get_communicator());

    AssertIndexRange(global_min, tria.n_global_levels());

    return global_min;
  }
} // namespace MGTools


// explicit instantiations
#include "mg_tools.inst"

DEAL_II_NAMESPACE_CLOSE