grid_tools_dof_handlers.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#include <deal.II/base/point.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/tensor.h>

#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/shared_tria.h>

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

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

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

#include <deal.II/hp/mapping_collection.h>
#include <deal.II/hp/dof_handler.h>

#include <deal.II/lac/full_matrix.h>

#include <algorithm>
#include <array>
#include <cmath>
#include <list>
#include <map>
#include <numeric>
#include <set>
#include <vector>


DEAL_II_NAMESPACE_OPEN

namespace GridTools
{
  template <int dim, template <int, int> class MeshType, int spacedim>
  unsigned int
  find_closest_vertex (const MeshType<dim,spacedim> &mesh,
                       const Point<spacedim>        &p,
                       const std::vector<bool>      &marked_vertices)
  {
    // first get the underlying
    // triangulation from the
    // mesh and determine vertices
    // and used vertices
    const Triangulation<dim, spacedim> &tria = mesh.get_triangulation();

    const std::vector< Point<spacedim> > &vertices = tria.get_vertices();

    Assert ( tria.get_vertices().size() == marked_vertices.size() || marked_vertices.size() ==0,
             ExcDimensionMismatch(tria.get_vertices().size(), marked_vertices.size()));

    // If p is an element of marked_vertices,
    // and q is that of used_Vertices,
    // the vector marked_vertices does NOT
    // contain unused vertices if p implies q.
    // I.e., if p is true q must be true
    // (if p is false, q could be false or true).
    // p implies q logic is encapsulated in ~p|q.
    Assert( marked_vertices.size()==0
            ||
            std::equal( marked_vertices.begin(),
                        marked_vertices.end(),
                        tria.get_used_vertices().begin(),
                        [](bool p, bool q)
    {
      return !p || q;
    }),
    ExcMessage("marked_vertices should be a subset of used vertices in the triangulation "
               "but marked_vertices contains one or more vertices that are not used vertices!") );

    // In addition, if a vector bools
    // is specified (marked_vertices)
    // marking all the vertices which
    // could be the potentially closest
    // vertex to the point, use it instead
    // of used vertices
    const std::vector<bool> &used     =
      (marked_vertices.size()==0) ? tria.get_used_vertices() : marked_vertices;

    // At the beginning, the first
    // used vertex is the closest one
    std::vector<bool>::const_iterator first =
      std::find(used.begin(), used.end(), true);

    // Assert that at least one vertex
    // is actually used
    Assert(first != used.end(), ExcInternalError());

    unsigned int best_vertex = std::distance(used.begin(), first);
    double       best_dist   = (p - vertices[best_vertex]).norm_square();

    // For all remaining vertices, test
    // whether they are any closer
    for (unsigned int j = best_vertex+1; j < vertices.size(); j++)
      if (used[j])
        {
          double dist = (p - vertices[j]).norm_square();
          if (dist < best_dist)
            {
              best_vertex = j;
              best_dist   = dist;
            }
        }

    return best_vertex;
  }



  template <int dim, template <int, int> class MeshType, int spacedim>
  unsigned int
  find_closest_vertex (const Mapping<dim,spacedim>  &mapping,
                       const MeshType<dim,spacedim> &mesh,
                       const Point<spacedim>        &p,
                       const std::vector<bool>      &marked_vertices)
  {
    // Take a shortcut in the simple case.
    if (mapping.preserves_vertex_locations() == true)
      return find_closest_vertex(mesh, p, marked_vertices);

    // first get the underlying
    // triangulation from the
    // mesh and determine vertices
    // and used vertices
    const Triangulation<dim, spacedim> &tria = mesh.get_triangulation();

    auto vertices = extract_used_vertices(tria, mapping);

    Assert ( tria.get_vertices().size() == marked_vertices.size() || marked_vertices.size() ==0,
             ExcDimensionMismatch(tria.get_vertices().size(), marked_vertices.size()));

    // If p is an element of marked_vertices,
    // and q is that of used_Vertices,
    // the vector marked_vertices does NOT
    // contain unused vertices if p implies q.
    // I.e., if p is true q must be true
    // (if p is false, q could be false or true).
    // p implies q logic is encapsulated in ~p|q.
    Assert( marked_vertices.size()==0
            ||
            std::equal( marked_vertices.begin(),
                        marked_vertices.end(),
                        tria.get_used_vertices().begin(),
                        [](bool p, bool q)
    {
      return !p || q;
    }),
    ExcMessage("marked_vertices should be a subset of used vertices in the triangulation "
               "but marked_vertices contains one or more vertices that are not used vertices!") );

    // Remove from the map unwanted elements.
    if (marked_vertices.size())
      for (auto it = vertices.begin(); it != vertices.end(); )
        {
          if (marked_vertices[it->first] == false)
            {
              vertices.erase(it++);
            }
          else
            {
              ++it;
            }
        }

    return find_closest_vertex(vertices, p);
  }



  template <int dim, template <int, int> class MeshType, int spacedim>
#ifndef _MSC_VER
  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
#else
  std::vector<typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type>
#endif
  find_cells_adjacent_to_vertex(const MeshType<dim,spacedim> &mesh,
                                const unsigned int            vertex)
  {
    // make sure that the given vertex is
    // an active vertex of the underlying
    // triangulation
    Assert(vertex < mesh.get_triangulation().n_vertices(),
           ExcIndexRange(0,mesh.get_triangulation().n_vertices(),vertex));
    Assert(mesh.get_triangulation().get_used_vertices()[vertex],
           ExcVertexNotUsed(vertex));

    // use a set instead of a vector
    // to ensure that cells are inserted only
    // once
    std::set<typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type> adjacent_cells;

    typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type
    cell = mesh.begin_active(),
    endc = mesh.end();

    // go through all active cells and look if the vertex is part of that cell
    //
    // in 1d, this is all we need to care about. in 2d/3d we also need to worry
    // that the vertex might be a hanging node on a face or edge of a cell; in
    // this case, we would want to add those cells as well on whose faces the
    // vertex is located but for which it is not a vertex itself.
    //
    // getting this right is a lot simpler in 2d than in 3d. in 2d, a hanging
    // node can only be in the middle of a face and we can query the neighboring
    // cell from the current cell. on the other hand, in 3d a hanging node
    // vertex can also be on an edge but there can be many other cells on
    // this edge and we can not access them from the cell we are currently
    // on.
    //
    // so, in the 3d case, if we run the algorithm as in 2d, we catch all
    // those cells for which the vertex we seek is on a *subface*, but we
    // miss the case of cells for which the vertex we seek is on a
    // sub-edge for which there is no corresponding sub-face (because the
    // immediate neighbor behind this face is not refined), see for example
    // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we
    // haven't yet found the vertex for the current cell we also need to
    // look at the mid-points of edges
    //
    // as a final note, deciding whether a neighbor is actually coarser is
    // simple in the case of isotropic refinement (we just need to look at
    // the level of the current and the neighboring cell). however, this
    // isn't so simple if we have used anisotropic refinement since then
    // the level of a cell is not indicative of whether it is coarser or
    // not than the current cell. ultimately, we want to add all cells on
    // which the vertex is, independent of whether they are coarser or
    // finer and so in the 2d case below we simply add *any* *active* neighbor.
    // in the worst case, we add cells multiple times to the adjacent_cells
    // list, but std::set throws out those cells already entered
    for (; cell != endc; ++cell)
      {
        for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; v++)
          if (cell->vertex_index(v) == vertex)
            {
              // OK, we found a cell that contains
              // the given vertex. We add it
              // to the list.
              adjacent_cells.insert(cell);

              // as explained above, in 2+d we need to check whether
              // this vertex is on a face behind which there is a
              // (possibly) coarser neighbor. if this is the case,
              // then we need to also add this neighbor
              if (dim >= 2)
                for (unsigned int vface = 0; vface < dim; vface++)
                  {
                    const unsigned int face =
                      GeometryInfo<dim>::vertex_to_face[v][vface];

                    if (!cell->at_boundary(face)
                        &&
                        cell->neighbor(face)->active())
                      {
                        // there is a (possibly) coarser cell behind a
                        // face to which the vertex belongs. the
                        // vertex we are looking at is then either a
                        // vertex of that coarser neighbor, or it is a
                        // hanging node on one of the faces of that
                        // cell. in either case, it is adjacent to the
                        // vertex, so add it to the list as well (if
                        // the cell was already in the list then the
                        // std::set makes sure that we get it only
                        // once)
                        adjacent_cells.insert (cell->neighbor(face));
                      }
                  }

              // in any case, we have found a cell, so go to the next cell
              goto next_cell;
            }

        // in 3d also loop over the edges
        if (dim >= 3)
          {
            for (unsigned int e=0; e<GeometryInfo<dim>::lines_per_cell; ++e)
              if (cell->line(e)->has_children())
                // the only place where this vertex could have been
                // hiding is on the mid-edge point of the edge we
                // are looking at
                if (cell->line(e)->child(0)->vertex_index(1) == vertex)
                  {
                    adjacent_cells.insert(cell);

                    // jump out of this tangle of nested loops
                    goto next_cell;
                  }
          }

        // in more than 3d we would probably have to do the same as
        // above also for even lower-dimensional objects
        Assert (dim <= 3, ExcNotImplemented());

        // move on to the next cell if we have found the
        // vertex on the current one
next_cell:
        ;
      }

    // if this was an active vertex then there needs to have been
    // at least one cell to which it is adjacent!
    Assert (adjacent_cells.size() > 0, ExcInternalError());

    // return the result as a vector, rather than the set we built above
    return
      std::vector<typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type>
      (adjacent_cells.begin(), adjacent_cells.end());
  }



  namespace
  {
    template <int dim, template <int, int> class MeshType, int spacedim>
    void find_active_cell_around_point_internal
    (const MeshType<dim,spacedim> &mesh,
#ifndef _MSC_VER
     std::set<typename MeshType<dim, spacedim>::active_cell_iterator> &searched_cells,
     std::set<typename MeshType<dim, spacedim>::active_cell_iterator> &adjacent_cells)
#else
     std::set<typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type> &searched_cells,
     std::set<typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type> &adjacent_cells)
#endif
    {
#ifndef _MSC_VER
      typedef typename MeshType<dim, spacedim>::active_cell_iterator cell_iterator;
#else
      typedef typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type cell_iterator;
#endif

      // update the searched cells
      searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end());
      // now we to collect all neighbors
      // of the cells in adjacent_cells we
      // have not yet searched.
      std::set<cell_iterator> adjacent_cells_new;

      typename std::set<cell_iterator>::const_iterator
      cell = adjacent_cells.begin(),
      endc = adjacent_cells.end();
      for (; cell != endc; ++cell)
        {
          std::vector<cell_iterator> active_neighbors;
          get_active_neighbors<MeshType<dim, spacedim> >(*cell, active_neighbors);
          for (unsigned int i=0; i<active_neighbors.size(); ++i)
            if (searched_cells.find(active_neighbors[i]) == searched_cells.end())
              adjacent_cells_new.insert(active_neighbors[i]);
        }
      adjacent_cells.clear();
      adjacent_cells.insert(adjacent_cells_new.begin(), adjacent_cells_new.end());
      if (adjacent_cells.size() == 0)
        {
          // we haven't found any other cell that would be a
          // neighbor of a previously found cell, but we know
          // that we haven't checked all cells yet. that means
          // that the domain is disconnected. in that case,
          // choose the first previously untouched cell we
          // can find
          cell_iterator it = mesh.begin_active();
          for ( ; it!=mesh.end(); ++it)
            if (searched_cells.find(it) == searched_cells.end())
              {
                adjacent_cells.insert(it);
                break;
              }
        }
    }
  }

  template <int dim, template <int, int> class MeshType, int spacedim>
#ifndef _MSC_VER
  typename MeshType<dim, spacedim>::active_cell_iterator
#else
  typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type
#endif
  find_active_cell_around_point (const MeshType<dim,spacedim> &mesh,
                                 const Point<spacedim>        &p,
                                 const std::vector<bool>      &marked_vertices)
  {
    return
      find_active_cell_around_point<dim,MeshType,spacedim>
      (StaticMappingQ1<dim,spacedim>::mapping,
       mesh, p, marked_vertices).first;
  }


  template <int dim, template <int, int> class MeshType, int spacedim>
#ifndef _MSC_VER
  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> >
#else
  std::pair<typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type, Point<dim> >
#endif
  find_active_cell_around_point (const Mapping<dim,spacedim>  &mapping,
                                 const MeshType<dim,spacedim> &mesh,
                                 const Point<spacedim>        &p,
                                 const std::vector<bool>      &marked_vertices)
  {
    typedef typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type active_cell_iterator;

    // The best distance is set to the
    // maximum allowable distance from
    // the unit cell; we assume a
    // max. deviation of 1e-10
    double best_distance = 1e-10;
    int    best_level = -1;
    std::pair<active_cell_iterator, Point<dim> > best_cell;

    // Find closest vertex and determine
    // all adjacent cells
    std::vector<active_cell_iterator> adjacent_cells_tmp
      = find_cells_adjacent_to_vertex(mesh,
                                      find_closest_vertex(mapping, mesh, p, marked_vertices));

    // Make sure that we have found
    // at least one cell adjacent to vertex.
    Assert(adjacent_cells_tmp.size()>0, ExcInternalError());

    // Copy all the cells into a std::set
    std::set<active_cell_iterator> adjacent_cells (adjacent_cells_tmp.begin(),
                                                   adjacent_cells_tmp.end());
    std::set<active_cell_iterator> searched_cells;

    // Determine the maximal number of cells
    // in the grid.
    // As long as we have not found
    // the cell and have not searched
    // every cell in the triangulation,
    // we keep on looking.
    const unsigned int n_active_cells = mesh.get_triangulation().n_active_cells();
    bool found = false;
    unsigned int cells_searched = 0;
    while (!found && cells_searched < n_active_cells)
      {
        typename std::set<active_cell_iterator>::const_iterator
        cell = adjacent_cells.begin(),
        endc = adjacent_cells.end();
        for (; cell != endc; ++cell)
          {
            try
              {
                const Point<dim> p_cell = mapping.transform_real_to_unit_cell(*cell, p);

                // calculate the infinity norm of
                // the distance vector to the unit cell.
                const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);

                // We compare if the point is inside the
                // unit cell (or at least not too far
                // outside). If it is, it is also checked
                // that the cell has a more refined state
                if ((dist < best_distance)
                    ||
                    ((dist == best_distance)
                     &&
                     ((*cell)->level() > best_level)))
                  {
                    found         = true;
                    best_distance = dist;
                    best_level    = (*cell)->level();
                    best_cell     = std::make_pair(*cell, p_cell);
                  }
              }
            catch (typename MappingQGeneric<dim,spacedim>::ExcTransformationFailed &)
              {
                // ok, the transformation
                // failed presumably
                // because the point we
                // are looking for lies
                // outside the current
                // cell. this means that
                // the current cell can't
                // be the cell around the
                // point, so just ignore
                // this cell and move on
                // to the next
              }
          }

        // update the number of cells searched
        cells_searched += adjacent_cells.size();

        // if the user provided a custom mask for vertices,
        // terminate the search without trying to expand the search
        // to all cells of the triangulation, as done below.
        if (marked_vertices.size() > 0)
          cells_searched = n_active_cells;

        // if we have not found the cell in
        // question and have not yet searched every
        // cell, we expand our search to
        // all the not already searched neighbors of
        // the cells in adjacent_cells. This is
        // what find_active_cell_around_point_internal
        // is for.
        if (!found && cells_searched < n_active_cells)
          {
            find_active_cell_around_point_internal<dim,MeshType,spacedim>
            (mesh, searched_cells, adjacent_cells);
          }
      }

    AssertThrow (best_cell.first.state() == IteratorState::valid,
                 ExcPointNotFound<spacedim>(p));

    return best_cell;
  }



  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  compute_active_cell_halo_layer
  (const MeshType                                                                    &mesh,
   const std::function<bool (const typename MeshType::active_cell_iterator &)> &predicate)
  {
    std::vector<typename MeshType::active_cell_iterator> active_halo_layer;
    std::vector<bool> locally_active_vertices_on_subdomain (mesh.get_triangulation().n_vertices(),
                                                            false);

    // Find the cells for which the predicate is true
    // These are the cells around which we wish to construct
    // the halo layer
    for (typename MeshType::active_cell_iterator
         cell = mesh.begin_active();
         cell != mesh.end(); ++cell)
      if (predicate(cell)) // True predicate --> Part of subdomain
        for (unsigned int v=0; v<GeometryInfo<MeshType::dimension>::vertices_per_cell; ++v)
          locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true;

    // Find the cells that do not conform to the predicate
    // but share a vertex with the selected subdomain
    // These comprise the halo layer
    for (typename MeshType::active_cell_iterator
         cell = mesh.begin_active();
         cell != mesh.end(); ++cell)
      if (!predicate(cell)) // False predicate --> Potential halo cell
        for (unsigned int v=0; v<GeometryInfo<MeshType::dimension>::vertices_per_cell; ++v)
          if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] == true)
            {
              active_halo_layer.push_back(cell);
              break;
            }

    return active_halo_layer;
  }



  template <class MeshType>
  std::vector<typename MeshType::cell_iterator>
  compute_cell_halo_layer_on_level
  (const MeshType                                                       &mesh,
   const std::function<bool (const typename MeshType::cell_iterator &)> &predicate,
   const unsigned int                                                    level)
  {
    std::vector<typename MeshType::cell_iterator> level_halo_layer;
    std::vector<bool> locally_active_vertices_on_level_subdomain (mesh.get_triangulation().n_vertices(),
        false);

    // Find the cells for which the predicate is true
    // These are the cells around which we wish to construct
    // the halo layer
    for (typename MeshType::cell_iterator
         cell = mesh.begin(level);
         cell != mesh.end(level); ++cell)
      if (predicate(cell)) // True predicate --> Part of subdomain
        for (unsigned int v=0; v<GeometryInfo<MeshType::dimension>::vertices_per_cell; ++v)
          locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] = true;

    // Find the cells that do not conform to the predicate
    // but share a vertex with the selected subdomain on that level
    // These comprise the halo layer
    for (typename MeshType::cell_iterator
         cell = mesh.begin(level);
         cell != mesh.end(level); ++cell)
      if (!predicate(cell)) // False predicate --> Potential halo cell
        for (unsigned int v=0; v<GeometryInfo<MeshType::dimension>::vertices_per_cell; ++v)
          if (locally_active_vertices_on_level_subdomain[cell->vertex_index(v)] == true)
            {
              level_halo_layer.push_back(cell);
              break;
            }

    return level_halo_layer;
  }


  namespace
  {
    template <class MeshType>
    bool
    contains_locally_owned_cells (const std::vector<typename MeshType::active_cell_iterator> &cells)
    {
      for (typename std::vector<typename MeshType::active_cell_iterator>::const_iterator
           it = cells.begin(); it != cells.end(); ++it)
        {
          if ((*it)->is_locally_owned())
            return true;
        }
      return false;
    }

    template <class MeshType>
    bool
    contains_artificial_cells (const std::vector<typename MeshType::active_cell_iterator> &cells)
    {
      for (typename std::vector<typename MeshType::active_cell_iterator>::const_iterator
           it = cells.begin(); it != cells.end(); ++it)
        {
          if ((*it)->is_artificial())
            return true;
        }
      return false;
    }
  }




  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  compute_ghost_cell_halo_layer (const MeshType &mesh)
  {
    std::function<bool (const typename MeshType::active_cell_iterator &)> predicate
      = IteratorFilters::LocallyOwnedCell();

    const std::vector<typename MeshType::active_cell_iterator>
    active_halo_layer = compute_active_cell_halo_layer (mesh, predicate);

    // Check that we never return locally owned or artificial cells
    // What is left should only be the ghost cells
    Assert(contains_locally_owned_cells<MeshType>(active_halo_layer) == false,
           ExcMessage("Halo layer contains locally owned cells"));
    Assert(contains_artificial_cells<MeshType>(active_halo_layer) == false,
           ExcMessage("Halo layer contains artificial cells"));

    return active_halo_layer;
  }



  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  compute_active_cell_layer_within_distance
  (const MeshType                                                                    &mesh,
   const std::function<bool (const typename MeshType::active_cell_iterator &)> &predicate,
   const double                                                                       layer_thickness)
  {
    std::vector<typename MeshType::active_cell_iterator> subdomain_boundary_cells, active_cell_layer_within_distance;
    std::vector<bool> vertices_outside_subdomain ( mesh.get_triangulation().n_vertices(),
                                                   false);

    const unsigned int spacedim = MeshType::space_dimension;

    unsigned int n_non_predicate_cells = 0; // Number of non predicate cells

    // Find the layer of cells for which predicate is true and that
    // are on the boundary with other cells. These are
    // subdomain boundary cells.

    // Find the cells for which the predicate is false
    // These are the cells which are around the predicate subdomain
    for ( typename MeshType::active_cell_iterator
          cell = mesh.begin_active();
          cell != mesh.end(); ++cell)
      if ( !predicate(cell)) // Negation of predicate --> Not Part of subdomain
        {
          for (unsigned int v=0; v<GeometryInfo<MeshType::dimension>::vertices_per_cell; ++v)
            vertices_outside_subdomain[cell->vertex_index(v)] = true;
          n_non_predicate_cells++;
        }

    // If all the active cells conform to the predicate
    // or if none of the active cells conform to the predicate
    // there is no active cell layer around the predicate
    // subdomain (within any distance)
    if ( n_non_predicate_cells == 0  || n_non_predicate_cells == mesh.get_triangulation().n_active_cells() )
      return std::vector<typename MeshType::active_cell_iterator>();

    // Find the cells that conform to the predicate
    // but share a vertex with the cell not in the predicate subdomain
    for ( typename MeshType::active_cell_iterator
          cell = mesh.begin_active();
          cell != mesh.end(); ++cell)
      if ( predicate(cell)) // True predicate --> Potential boundary cell of the subdomain
        for (unsigned int v=0; v<GeometryInfo<MeshType::dimension>::vertices_per_cell; ++v)
          if (vertices_outside_subdomain[cell->vertex_index(v)] == true)
            {
              subdomain_boundary_cells.push_back(cell);
              break; // No need to go through remaining vertices
            }

    // To cheaply filter out some cells located far away from the predicate subdomain,
    // get the bounding box of the predicate subdomain.
    std::pair< Point<spacedim>, Point<spacedim> > bounding_box = compute_bounding_box( mesh,
        predicate );

    // DOUBLE_EPSILON to compare really close double values
    const double &DOUBLE_EPSILON = 100.*std::numeric_limits<double>::epsilon();

    // Add layer_thickness to the bounding box
    for ( unsigned int d=0; d<spacedim; ++d)
      {
        bounding_box.first[d]  -= (layer_thickness+DOUBLE_EPSILON); // minp
        bounding_box.second[d] += (layer_thickness+DOUBLE_EPSILON); // maxp
      }

    std::vector<Point<spacedim> > subdomain_boundary_cells_centers; // cache all the subdomain boundary cells centers here
    std::vector<double> subdomain_boundary_cells_radii; // cache all the subdomain boundary cells radii
    subdomain_boundary_cells_centers.reserve (subdomain_boundary_cells.size());
    subdomain_boundary_cells_radii.reserve (subdomain_boundary_cells.size());
    // compute cell radius for each boundary cell of the predicate subdomain
    for ( typename std::vector<typename MeshType::active_cell_iterator>::const_iterator
          subdomain_boundary_cell_iterator  = subdomain_boundary_cells.begin();
          subdomain_boundary_cell_iterator != subdomain_boundary_cells.end(); ++subdomain_boundary_cell_iterator )
      {
        const std::pair<Point<spacedim>, double> &
        subdomain_boundary_cell_enclosing_ball = (*subdomain_boundary_cell_iterator)->enclosing_ball();

        subdomain_boundary_cells_centers.push_back( subdomain_boundary_cell_enclosing_ball.first);
        subdomain_boundary_cells_radii.push_back( subdomain_boundary_cell_enclosing_ball.second);
      }
    AssertThrow( subdomain_boundary_cells_radii.size() == subdomain_boundary_cells_centers.size(),
                 ExcInternalError());

    // Find the cells that are within layer_thickness of predicate subdomain boundary
    // distance but are inside the extended bounding box.
    // Most cells might be outside the extended bounding box, so we could skip them.
    // Those cells that are inside the extended bounding box but are not part of the
    // predicate subdomain are possible candidates to be within the distance to the
    // boundary cells of the predicate subdomain.
    for ( typename MeshType::active_cell_iterator
          cell = mesh.begin_active();
          cell != mesh.end(); ++cell)
      {
        // Ignore all the cells that are in the predicate subdomain
        if ( predicate(cell))
          continue;

        const std::pair<Point<spacedim>, double> &cell_enclosing_ball
          = cell->enclosing_ball();

        const Point<spacedim> &cell_enclosing_ball_center = cell_enclosing_ball.first;
        const double &cell_enclosing_ball_radius = cell_enclosing_ball.second;

        bool cell_inside = true; // reset for each cell

        for (unsigned int d = 0; d < spacedim; ++d)
          cell_inside &= (cell_enclosing_ball_center[d] + cell_enclosing_ball_radius > bounding_box.first[d])
                         && (cell_enclosing_ball_center[d] - cell_enclosing_ball_radius < bounding_box.second[d]);
        // cell_inside is true if its enclosing ball intersects the extended bounding box

        // Ignore all the cells that are outside the extended bounding box
        if (cell_inside)
          for (unsigned int i =0; i< subdomain_boundary_cells_radii.size(); ++i)
            if ( cell_enclosing_ball_center.distance_square(subdomain_boundary_cells_centers[i])
                 <  Utilities::fixed_power<2>( cell_enclosing_ball_radius +
                                               subdomain_boundary_cells_radii[i] +
                                               layer_thickness + DOUBLE_EPSILON ))
              {
                active_cell_layer_within_distance.push_back(cell);
                break; // Exit the loop checking all the remaining subdomain boundary cells
              }

      }
    return active_cell_layer_within_distance;
  }



  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  compute_ghost_cell_layer_within_distance ( const MeshType &mesh, const double layer_thickness)
  {
    IteratorFilters::LocallyOwnedCell locally_owned_cell_predicate;
    std::function<bool (const typename MeshType::active_cell_iterator &)> predicate (locally_owned_cell_predicate);

    const std::vector<typename MeshType::active_cell_iterator>
    ghost_cell_layer_within_distance = compute_active_cell_layer_within_distance (mesh, predicate, layer_thickness);

    // Check that we never return locally owned or artificial cells
    // What is left should only be the ghost cells
    Assert(contains_locally_owned_cells<MeshType>(ghost_cell_layer_within_distance) == false,
           ExcMessage("Ghost cells within layer_thickness contains locally owned cells."));
    Assert(contains_artificial_cells<MeshType>(ghost_cell_layer_within_distance) == false,
           ExcMessage("Ghost cells within layer_thickness contains artificial cells."
                      "The function compute_ghost_cell_layer_within_distance "
                      "is probably called while using parallel::distributed::Triangulation. "
                      "In such case please refer to the description of this function."));

    return ghost_cell_layer_within_distance;
  }



  template < class MeshType>
  std::pair< Point<MeshType::space_dimension>, Point<MeshType::space_dimension> >
  compute_bounding_box
  ( const MeshType                                                                    &mesh,
    const std::function<bool (const typename MeshType::active_cell_iterator &)> &predicate )
  {
    std::vector<bool> locally_active_vertices_on_subdomain (mesh.get_triangulation().n_vertices(),
                                                            false);

    const unsigned int spacedim = MeshType::space_dimension;

    // Two extreme points can define the bounding box
    // around the active cells that conform to the given predicate.
    Point<MeshType::space_dimension> maxp, minp;

    // initialize minp and maxp with the first predicate cell center
    for ( typename MeshType::active_cell_iterator
          cell = mesh.begin_active();
          cell != mesh.end(); ++cell)
      if ( predicate(cell))
        {
          minp = cell->center();
          maxp = cell->center();
          break;
        }

    // Run through all the cells to check if it belongs to predicate domain,
    // if it belongs to the predicate domain, extend the bounding box.
    for ( typename MeshType::active_cell_iterator
          cell = mesh.begin_active();
          cell != mesh.end(); ++cell)
      if (predicate(cell)) // True predicate --> Part of subdomain
        for (unsigned int v=0; v<GeometryInfo<MeshType::dimension>::vertices_per_cell; ++v)
          if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] == false)
            {
              locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true;
              for ( unsigned int d=0; d<spacedim; ++d)
                {
                  minp[d] = std::min( minp[d], cell->vertex(v)[d]);
                  maxp[d] = std::max( maxp[d], cell->vertex(v)[d]);
                }
            }

    return std::make_pair( minp, maxp );
  }



  template <typename MeshType>
  std::list<std::pair<typename MeshType::cell_iterator,
      typename MeshType::cell_iterator> >
      get_finest_common_cells (const MeshType &mesh_1,
                               const MeshType &mesh_2)
  {
    Assert (have_same_coarse_mesh (mesh_1, mesh_2),
            ExcMessage ("The two meshes must be represent triangulations that "
                        "have the same coarse meshes"));

    // the algorithm goes as follows:
    // first, we fill a list with pairs
    // of iterators common to the two
    // meshes on the coarsest
    // level. then we traverse the
    // list; each time, we find a pair
    // of iterators for which both
    // correspond to non-active cells,
    // we delete this item and push the
    // pairs of iterators to their
    // children to the back. if these
    // again both correspond to
    // non-active cells, we will get to
    // the later on for further
    // consideration
    typedef
    std::list<std::pair<typename MeshType::cell_iterator,
        typename MeshType::cell_iterator> >
        CellList;

    CellList cell_list;

    // first push the coarse level cells
    typename MeshType::cell_iterator
    cell_1 = mesh_1.begin(0),
    cell_2 = mesh_2.begin(0);
    for (; cell_1 != mesh_1.end(0); ++cell_1, ++cell_2)
      cell_list.emplace_back (cell_1, cell_2);

    // then traverse list as described
    // above
    typename CellList::iterator cell_pair = cell_list.begin();
    while (cell_pair != cell_list.end())
      {
        // if both cells in this pair
        // have children, then erase
        // this element and push their
        // children instead
        if (cell_pair->first->has_children()
            &&
            cell_pair->second->has_children())
          {
            Assert(cell_pair->first->refinement_case()==
                   cell_pair->second->refinement_case(), ExcNotImplemented());
            for (unsigned int c=0; c<cell_pair->first->n_children(); ++c)
              cell_list.emplace_back (cell_pair->first->child(c),
                                      cell_pair->second->child(c));

            // erasing an iterator
            // keeps other iterators
            // valid, so already
            // advance the present
            // iterator by one and then
            // delete the element we've
            // visited before
            const typename CellList::iterator previous_cell_pair = cell_pair;
            ++cell_pair;

            cell_list.erase (previous_cell_pair);
          }
        else
          // both cells are active, do
          // nothing
          ++cell_pair;
      }

    // just to make sure everything is ok,
    // validate that all pairs have at least one
    // active iterator or have different
    // refinement_cases
    for (cell_pair = cell_list.begin(); cell_pair != cell_list.end(); ++cell_pair)
      Assert (cell_pair->first->active()
              ||
              cell_pair->second->active()
              ||
              (cell_pair->first->refinement_case()
               != cell_pair->second->refinement_case()),
              ExcInternalError());

    return cell_list;
  }



  template <int dim, int spacedim>
  bool
  have_same_coarse_mesh (const Triangulation<dim, spacedim> &mesh_1,
                         const Triangulation<dim, spacedim> &mesh_2)
  {
    // make sure the two meshes have
    // the same number of coarse cells
    if (mesh_1.n_cells (0) != mesh_2.n_cells (0))
      return false;

    // if so, also make sure they have
    // the same vertices on the cells
    // of the coarse mesh
    typename Triangulation<dim, spacedim>::cell_iterator
    cell_1 = mesh_1.begin(0),
    cell_2 = mesh_2.begin(0),
    endc   = mesh_1.end(0);
    for (; cell_1!=endc; ++cell_1, ++cell_2)
      for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
        if (cell_1->vertex(v) != cell_2->vertex(v))
          return false;

    // if we've gotten through all
    // this, then the meshes really
    // seem to have a common coarse
    // mesh
    return true;
  }



  template <typename MeshType>
  bool
  have_same_coarse_mesh (const MeshType &mesh_1,
                         const MeshType &mesh_2)
  {
    return have_same_coarse_mesh (mesh_1.get_triangulation(),
                                  mesh_2.get_triangulation());
  }



  template <int dim, int spacedim>
  std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<dim> >
  find_active_cell_around_point (const hp::MappingCollection<dim,spacedim>   &mapping,
                                 const hp::DoFHandler<dim,spacedim> &mesh,
                                 const Point<spacedim>     &p)
  {
    Assert ((mapping.size() == 1) ||
            (mapping.size() == mesh.get_fe_collection().size()),
            ExcMessage ("Mapping collection needs to have either size 1 "
                        "or size equal to the number of elements in "
                        "the FECollection."));

    typedef typename hp::DoFHandler<dim,spacedim>::active_cell_iterator cell_iterator;

    std::pair<cell_iterator, Point<dim> > best_cell;
    //If we have only one element in the MappingCollection,
    //we use find_active_cell_around_point using only one
    //mapping.
    if (mapping.size() == 1)
      best_cell = find_active_cell_around_point(mapping[0], mesh, p);
    else
      {


        // The best distance is set to the
        // maximum allowable distance from
        // the unit cell; we assume a
        // max. deviation of 1e-10
        double best_distance = 1e-10;
        int    best_level = -1;


        // Find closest vertex and determine
        // all adjacent cells
        unsigned int vertex = find_closest_vertex(mesh, p);

        std::vector<cell_iterator> adjacent_cells_tmp =
          find_cells_adjacent_to_vertex(mesh, vertex);

        // Make sure that we have found
        // at least one cell adjacent to vertex.
        Assert(adjacent_cells_tmp.size()>0, ExcInternalError());

        // Copy all the cells into a std::set
        std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end());
        std::set<cell_iterator> searched_cells;

        // Determine the maximal number of cells
        // in the grid.
        // As long as we have not found
        // the cell and have not searched
        // every cell in the triangulation,
        // we keep on looking.
        const unsigned int n_cells = mesh.get_triangulation().n_cells();
        bool found = false;
        unsigned int cells_searched = 0;
        while (!found && cells_searched < n_cells)
          {
            typename std::set<cell_iterator>::const_iterator
            cell = adjacent_cells.begin(),
            endc = adjacent_cells.end();
            for (; cell != endc; ++cell)
              {
                try
                  {
                    const Point<dim> p_cell = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p);


                    // calculate the infinity norm of
                    // the distance vector to the unit cell.
                    const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);

                    // We compare if the point is inside the
                    // unit cell (or at least not too far
                    // outside). If it is, it is also checked
                    // that the cell has a more refined state
                    if (dist < best_distance ||
                        (dist == best_distance && (*cell)->level() > best_level))
                      {
                        found       = true;
                        best_distance = dist;
                        best_level    = (*cell)->level();
                        best_cell     = std::make_pair(*cell, p_cell);
                      }
                  }
                catch (typename MappingQGeneric<dim,spacedim>::ExcTransformationFailed &)
                  {
                    // ok, the transformation
                    // failed presumably
                    // because the point we
                    // are looking for lies
                    // outside the current
                    // cell. this means that
                    // the current cell can't
                    // be the cell around the
                    // point, so just ignore
                    // this cell and move on
                    // to the next
                  }
              }
            //update the number of cells searched
            cells_searched += adjacent_cells.size();
            // if we have not found the cell in
            // question and have not yet searched every
            // cell, we expand our search to
            // all the not already searched neighbors of
            // the cells in adjacent_cells.
            if (!found && cells_searched < n_cells)
              {
                find_active_cell_around_point_internal<dim,hp::DoFHandler,spacedim>
                (mesh, searched_cells, adjacent_cells);
              }

          }
      }

    AssertThrow (best_cell.first.state() == IteratorState::valid,
                 ExcPointNotFound<spacedim>(p));

    return best_cell;
  }


  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell)
  {
    Assert (cell->is_locally_owned(),
            ExcMessage ("This function only makes sense if the cell for "
                        "which you are asking for a patch, is locally "
                        "owned."));

    std::vector<typename MeshType::active_cell_iterator> patch;
    patch.push_back (cell);
    for (unsigned int face_number=0; face_number<GeometryInfo<MeshType::dimension>::faces_per_cell; ++face_number)
      if (cell->face(face_number)->at_boundary()==false)
        {
          if (cell->neighbor(face_number)->has_children() == false)
            patch.push_back (cell->neighbor(face_number));
          else
            // the neighbor is refined. in 2d/3d, we can simply ask for the children
            // of the neighbor because they can not be further refined and,
            // consequently, the children is active
            if (MeshType::dimension > 1)
              {
                for (unsigned int subface=0; subface<cell->face(face_number)->n_children(); ++subface)
                  patch.push_back (cell->neighbor_child_on_subface (face_number, subface));
              }
            else
              {
                // in 1d, we need to work a bit harder: iterate until we find
                // the child by going from cell to child to child etc
                typename MeshType::cell_iterator neighbor
                  = cell->neighbor (face_number);
                while (neighbor->has_children())
                  neighbor = neighbor->child(1-face_number);

                Assert (neighbor->neighbor(1-face_number) == cell, ExcInternalError());
                patch.push_back (neighbor);
              }
        }
    return patch;
  }



  template <class Container>
  std::vector<typename Container::cell_iterator>
  get_cells_at_coarsest_common_level (
    const std::vector<typename Container::active_cell_iterator> &patch)
  {
    Assert (patch.size() > 0, ExcMessage("Vector containing patch cells should not be an empty vector!"));
    // In order to extract the set of cells with the coarsest common level from the give vector of cells:
    // First it finds the number associated with the minimum level of refinmenet, namely "min_level"
    int min_level = patch[0]->level();

    for (unsigned int i=0; i<patch.size(); ++i)
      min_level = std::min (min_level, patch[i]->level() );
    std::set<typename Container::cell_iterator>  uniform_cells;
    typename std::vector<typename Container::active_cell_iterator>::const_iterator  patch_cell;
    // it loops through all cells of the input vector
    for (patch_cell=patch.begin(); patch_cell!=patch.end () ; ++patch_cell)
      {
        // If the refinement level of each cell i the loop be equal to the min_level, so that
        // that cell inserted into the set of uniform_cells, as the set of cells with the coarsest common refinement level
        if ((*patch_cell)->level() == min_level)
          uniform_cells.insert (*patch_cell);
        else
          // If not, it asks for the parent of the cell, until it finds the parent cell
          // with the refinement level equal to the min_level and inserts that parent cell into the
          // the set of uniform_cells, as the set of cells with the coarsest common refinement level.
          {
            typename Container::cell_iterator parent = *patch_cell;

            while (parent->level() > min_level)
              parent = parent-> parent();
            uniform_cells.insert (parent);
          }
      }

    return std::vector<typename Container::cell_iterator> (uniform_cells.begin(),
                                                           uniform_cells.end());
  }



  template <class Container>
  void build_triangulation_from_patch(const std::vector<typename Container::active_cell_iterator> &patch,
                                      Triangulation<Container::dimension,Container::space_dimension> &local_triangulation,
                                      std::map<typename Triangulation<Container::dimension,Container::space_dimension>::active_cell_iterator,
                                      typename Container::active_cell_iterator> &patch_to_global_tria_map)

  {
    const std::vector<typename Container::cell_iterator> uniform_cells =
      get_cells_at_coarsest_common_level <Container> (patch);
    // First it creates triangulation from the vector of "uniform_cells"
    local_triangulation.clear();
    std::vector<Point<Container::space_dimension> > vertices;
    const unsigned int n_uniform_cells=uniform_cells.size();
    std::vector<CellData<Container::dimension> > cells(n_uniform_cells);
    unsigned int k=0;// for enumerating cells
    unsigned int i=0;// for enumerating vertices
    typename std::vector<typename Container::cell_iterator>::const_iterator uniform_cell;
    for (uniform_cell=uniform_cells.begin(); uniform_cell!=uniform_cells.end(); ++uniform_cell)
      {
        for (unsigned int v=0; v<GeometryInfo<Container::dimension>::vertices_per_cell; ++v)
          {
            Point<Container::space_dimension> position=(*uniform_cell)->vertex (v);
            bool repeat_vertex=false;

            for (unsigned int m=0; m<i; ++m)
              {
                if (position == vertices[m])
                  {
                    repeat_vertex=true;
                    cells[k].vertices[v]=m;
                    break;
                  }
              }
            if (repeat_vertex==false)
              {
                vertices.push_back(position);
                cells[k].vertices[v]=i;
                i=i+1;
              }

          }//for vertices_per_cell
        k=k+1;
      }
    local_triangulation.create_triangulation(vertices,cells,SubCellData());
    Assert (local_triangulation.n_active_cells() == uniform_cells.size(), ExcInternalError());
    local_triangulation.clear_user_flags ();
    unsigned int index=0;
    // Create a map between cells of class DoFHandler into class Triangulation
    std::map<typename Triangulation<Container::dimension,Container::space_dimension>::cell_iterator,
        typename Container::cell_iterator> patch_to_global_tria_map_tmp;
    for (typename Triangulation<Container::dimension,Container::space_dimension>::cell_iterator coarse_cell = local_triangulation.begin();
         coarse_cell != local_triangulation.end(); ++coarse_cell, ++index)
      {
        patch_to_global_tria_map_tmp.insert (std::make_pair(coarse_cell, uniform_cells[index]));
        // To ensure that the cells with the same coordinates (here, we compare their centers) are mapped into each other.

        Assert(coarse_cell->center().distance( uniform_cells[index]->center())<=1e-15*coarse_cell->diameter(),
               ExcInternalError());
      }
    bool refinement_necessary;
    // In this loop we start to do refinement on the above coarse triangulation to reach
    // to the same level of refinement as the patch cells are really on
    do
      {
        refinement_necessary = false;
        for (typename Triangulation<Container::dimension,Container::space_dimension>::active_cell_iterator
             active_tria_cell = local_triangulation.begin_active();
             active_tria_cell != local_triangulation.end(); ++active_tria_cell)
          {
            if (patch_to_global_tria_map_tmp[active_tria_cell]->has_children())
              {
                active_tria_cell -> set_refine_flag();
                refinement_necessary = true;
              }
            else for (unsigned int i=0; i<patch.size(); ++i)
                {
                  // Even though vertices may not be exactly the same, the
                  // appropriate cells will match since == for TriAccessors
                  // checks only cell level and index.
                  if (patch_to_global_tria_map_tmp[active_tria_cell]==patch[i])
                    {
                      // adjust the cell vertices of the local_triangulation to
                      // match cell vertices of the global triangulation
                      for (unsigned int v=0; v<GeometryInfo<Container::dimension>::vertices_per_cell; ++v)
                        active_tria_cell->vertex(v) = patch[i]->vertex(v);

                      Assert(active_tria_cell->center().distance(patch_to_global_tria_map_tmp[active_tria_cell]->center())
                             <=1e-15*active_tria_cell->diameter(), ExcInternalError());

                      active_tria_cell->set_user_flag();
                      break;
                    }
                }
          }

        if (refinement_necessary)
          {
            local_triangulation.execute_coarsening_and_refinement ();

            for (typename Triangulation<Container::dimension,Container::space_dimension>::cell_iterator
                 cell = local_triangulation.begin();
                 cell != local_triangulation.end(); ++cell)
              {

                if (patch_to_global_tria_map_tmp.find(cell)!=patch_to_global_tria_map_tmp.end())
                  {
                    if (cell-> has_children())
                      {
                        // Note: Since the cell got children, then it should not be in the map anymore
                        // children may be added into the map, instead

                        // these children may not yet be in the map
                        for (unsigned int c=0; c<cell->n_children(); ++c)
                          {
                            if (patch_to_global_tria_map_tmp.find(cell->child(c)) ==
                                patch_to_global_tria_map_tmp.end())
                              {
                                patch_to_global_tria_map_tmp.insert (std::make_pair(cell->child(c),
                                                                                    patch_to_global_tria_map_tmp[cell]->child(c)));

                                // One might be tempted to assert that the cell
                                // being added here has the same center as the
                                // equivalent cell in the global triangulation,
                                // but it may not be the case.  For triangulations
                                // that have been perturbed or smoothed, the cell
                                // indices and levels may be the same, but the
                                // vertex locations may not.  We adjust
                                // the vertices of the cells that have no
                                // children (ie the active cells) to be
                                // consistent with the global triangulation
                                // later on and add assertions at that time
                                // to guarantee the cells in the
                                // local_triangulation are physically at the same
                                // locations of the cells in the patch of the
                                // global triangulation.

                              }
                          }
                        // The parent cell whose children were added
                        // into the map should be deleted from the map
                        patch_to_global_tria_map_tmp.erase(cell);
                      }
                  }
              }
          }

      }
    while (refinement_necessary);


    // Last assertion check to make sure we have the right cells and centers
    // in the map, and hence the correct vertices of the triangulation
    for (typename Triangulation<Container::dimension,Container::space_dimension>::cell_iterator
         cell = local_triangulation.begin();
         cell != local_triangulation.end(); ++cell)
      {
        if (cell->user_flag_set() )
          {
            Assert(patch_to_global_tria_map_tmp.find(cell) != patch_to_global_tria_map_tmp.end(),
                   ExcInternalError() );

            Assert(cell->center().distance( patch_to_global_tria_map_tmp[cell]->center())<=1e-15*cell->diameter(),
                   ExcInternalError());
          }
      }


    typename std::map<typename Triangulation<Container::dimension,Container::space_dimension>::cell_iterator,
             typename Container::cell_iterator>::iterator map_tmp_it =
               patch_to_global_tria_map_tmp.begin(),map_tmp_end = patch_to_global_tria_map_tmp.end();
    // Now we just need to take the temporary map of pairs of type cell_iterator "patch_to_global_tria_map_tmp"
    // making pair of active_cell_iterators so that filling out the final map "patch_to_global_tria_map"
    for (; map_tmp_it!=map_tmp_end; ++map_tmp_it)
      patch_to_global_tria_map[map_tmp_it->first] = map_tmp_it->second;
  }




  template <class DoFHandlerType>
  std::map< types::global_dof_index,std::vector<typename DoFHandlerType::active_cell_iterator> >
  get_dof_to_support_patch_map(DoFHandlerType &dof_handler)
  {

    // This is the map from global_dof_index to
    // a set of cells on patch.  We first map into
    // a set because it is very likely that we
    // will attempt to add a cell more than once
    // to a particular patch and we want to preserve
    // uniqueness of cell iterators. std::set does this
    // automatically for us.  Later after it is all
    // constructed, we will copy to a map of vectors
    // since that is the preferred output for other
    // functions.
    std::map< types::global_dof_index,std::set<typename DoFHandlerType::active_cell_iterator> > dof_to_set_of_cells_map;

    std::vector<types::global_dof_index> local_dof_indices;
    std::vector<types::global_dof_index> local_face_dof_indices;
    std::vector<types::global_dof_index> local_line_dof_indices;

    // a place to save the dof_handler user flags and restore them later
    // to maintain const of dof_handler.
    std::vector<bool> user_flags;


    // in 3d, we need pointers from active lines to the
    // active parent lines, so we construct it as needed.
    std::map<typename DoFHandlerType::active_line_iterator, typename DoFHandlerType::line_iterator > lines_to_parent_lines_map;
    if (DoFHandlerType::dimension == 3)
      {

        // save user flags as they will be modified and then later restored
        dof_handler.get_triangulation().save_user_flags(user_flags);
        const_cast<::Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &>(dof_handler.get_triangulation()).clear_user_flags ();


        typename DoFHandlerType::active_cell_iterator cell = dof_handler.begin_active(),
                                                      endc = dof_handler.end();
        for (; cell!=endc; ++cell)
          {
            // We only want lines that are locally_relevant
            // although it doesn't hurt to have lines that
            // are children of ghost cells since there are
            // few and we don't have to use them.
            if (cell->is_artificial() == false)
              {
                for (unsigned int l=0; l<GeometryInfo<DoFHandlerType::dimension>::lines_per_cell; ++l)
                  if (cell->line(l)->has_children())
                    for (unsigned int c=0; c<cell->line(l)->n_children(); ++c)
                      {
                        lines_to_parent_lines_map[cell->line(l)->child(c)] = cell->line(l);
                        // set flags to know that child
                        // line has an active parent.
                        cell->line(l)->child(c)->set_user_flag();
                      }
              }
          }
      }


    // We loop through all cells and add cell to the
    // map for the dofs that it immediately touches
    // and then account for all the other dofs of
    // which it is a part, mainly the ones that must
    // be added on account of adaptivity hanging node
    // constraints.
    typename DoFHandlerType::active_cell_iterator cell = dof_handler.begin_active(),
                                                  endc = dof_handler.end();
    for (; cell!=endc; ++cell)
      {
        // Need to loop through all cells that could
        // be in the patch of dofs on locally_owned
        // cells including ghost cells
        if (cell->is_artificial() == false)
          {
            const unsigned int n_dofs_per_cell = cell->get_fe().dofs_per_cell;
            local_dof_indices.resize(n_dofs_per_cell);

            // Take care of adding cell pointer to each
            // dofs that exists on cell.
            cell->get_dof_indices(local_dof_indices);
            for (unsigned int i=0; i< n_dofs_per_cell; ++i )
              dof_to_set_of_cells_map[local_dof_indices[i]].insert(cell);

            // In the case of the adjacent cell (over
            // faces or edges) being more refined, we
            // want to add all of the children to the
            // patch since the support function at that
            // dof could be non-zero along that entire
            // face (or line).

            // Take care of dofs on neighbor faces
            for (unsigned int f=0; f<GeometryInfo<DoFHandlerType::dimension>::faces_per_cell; ++f)
              {
                if (cell->face(f)->has_children())
                  {
                    for (unsigned int c=0; c<cell->face(f)->n_children(); ++c)
                      {
                        //  Add cell to dofs of all subfaces
                        //
                        //   *-------------------*----------*---------*
                        //   |                   | add cell |         |
                        //   |                   |<- to dofs|         |
                        //   |                   |of subface|         |
                        //   |        cell       *----------*---------*
                        //   |                   | add cell |         |
                        //   |                   |<- to dofs|         |
                        //   |                   |of subface|         |
                        //   *-------------------*----------*---------*
                        //
                        Assert (cell->face(f)->child(c)->has_children() == false, ExcInternalError());

                        const unsigned int n_dofs_per_face = cell->get_fe().dofs_per_face;
                        local_face_dof_indices.resize(n_dofs_per_face);

                        cell->face(f)->child(c)->get_dof_indices(local_face_dof_indices);
                        for (unsigned int i=0; i< n_dofs_per_face; ++i )
                          dof_to_set_of_cells_map[local_face_dof_indices[i]].insert(cell);
                      }
                  }
                else if ((cell->face(f)->at_boundary() == false) && (cell->neighbor_is_coarser(f)))
                  {

                    // Add cell to dofs of parent face and all
                    // child faces of parent face
                    //
                    //   *-------------------*----------*---------*
                    //   |                   |          |         |
                    //   |                   |   cell   |         |
                    //   |      add cell     |          |         |
                    //   |      to dofs   -> *----------*---------*
                    //   |      of parent    | add cell |         |
                    //   |       face        |<- to dofs|         |
                    //   |                   |of subface|         |
                    //   *-------------------*----------*---------*
                    //

                    // Add cell to all dofs of parent face
                    std::pair<unsigned int, unsigned int> neighbor_face_no_subface_no = cell->neighbor_of_coarser_neighbor(f);
                    unsigned int face_no = neighbor_face_no_subface_no.first;
                    unsigned int subface = neighbor_face_no_subface_no.second;

                    const unsigned int n_dofs_per_face = cell->get_fe().dofs_per_face;
                    local_face_dof_indices.resize(n_dofs_per_face);

                    cell->neighbor(f)->face(face_no)->get_dof_indices(local_face_dof_indices);
                    for (unsigned int i=0; i< n_dofs_per_face; ++i )
                      dof_to_set_of_cells_map[local_face_dof_indices[i]].insert(cell);

                    // Add cell to all dofs of children of
                    // parent face
                    for (unsigned int c=0; c<cell->neighbor(f)->face(face_no)->n_children(); ++c)
                      {
                        if (c != subface) // don't repeat work on dofs of original cell
                          {
                            const unsigned int n_dofs_per_face = cell->get_fe().dofs_per_face;
                            local_face_dof_indices.resize(n_dofs_per_face);

                            Assert (cell->neighbor(f)->face(face_no)->child(c)->has_children() == false, ExcInternalError());
                            cell->neighbor(f)->face(face_no)->child(c)->get_dof_indices(local_face_dof_indices);
                            for (unsigned int i=0; i<n_dofs_per_face; ++i )
                              dof_to_set_of_cells_map[local_face_dof_indices[i]].insert(cell);
                          }
                      }
                  }
              }


            // If 3d, take care of dofs on lines in the
            // same pattern as faces above. That is, if
            // a cell's line has children, distribute
            // cell to dofs of children of line,  and
            // if cell's line has an active parent, then
            // distribute cell to dofs on parent line
            // and dofs on all children of parent line.
            if (DoFHandlerType::dimension == 3)
              {
                for (unsigned int l=0; l<GeometryInfo<DoFHandlerType::dimension>::lines_per_cell; ++l)
                  {
                    if (cell->line(l)->has_children())
                      {
                        for (unsigned int c=0; c<cell->line(l)->n_children(); ++c)
                          {
                            Assert (cell->line(l)->child(c)->has_children() == false, ExcInternalError());

                            // dofs_per_line returns number of dofs
                            // on line not including the vertices of the line.
                            const unsigned int n_dofs_per_line = 2*cell->get_fe().dofs_per_vertex
                                                                 + cell->get_fe().dofs_per_line;
                            local_line_dof_indices.resize(n_dofs_per_line);

                            cell->line(l)->child(c)->get_dof_indices(local_line_dof_indices);
                            for (unsigned int i=0; i<n_dofs_per_line; ++i )
                              dof_to_set_of_cells_map[local_line_dof_indices[i]].insert(cell);
                          }
                      }
                    // user flag was set above to denote that
                    // an active parent line exists so add
                    // cell to dofs of parent and all it's
                    // children
                    else if (cell->line(l)->user_flag_set() == true)
                      {
                        typename DoFHandlerType::line_iterator parent_line = lines_to_parent_lines_map[cell->line(l)];
                        Assert (parent_line->has_children(), ExcInternalError() );

                        // dofs_per_line returns number of dofs
                        // on line not including the vertices of the line.
                        const unsigned int n_dofs_per_line = 2*cell->get_fe().dofs_per_vertex
                                                             + cell->get_fe().dofs_per_line;
                        local_line_dof_indices.resize(n_dofs_per_line);

                        parent_line->get_dof_indices(local_line_dof_indices);
                        for (unsigned int i=0; i<n_dofs_per_line; ++i )
                          dof_to_set_of_cells_map[local_line_dof_indices[i]].insert(cell);

                        for (unsigned int c=0; c<parent_line->n_children(); ++c)
                          {
                            Assert (parent_line->child(c)->has_children() == false, ExcInternalError());

                            const unsigned int n_dofs_per_line = 2*cell->get_fe().dofs_per_vertex
                                                                 + cell->get_fe().dofs_per_line;
                            local_line_dof_indices.resize(n_dofs_per_line);

                            parent_line->child(c)->get_dof_indices(local_line_dof_indices);
                            for (unsigned int i=0; i<n_dofs_per_line; ++i )
                              dof_to_set_of_cells_map[local_line_dof_indices[i]].insert(cell);
                          }


                      }
                  } // for lines l
              }// if DoFHandlerType::dimension == 3
          }// if cell->is_locally_owned()
      }// for cells


    if (DoFHandlerType::dimension == 3)
      {
        // finally, restore user flags that were changed above
        // to when we constructed the pointers to parent of lines
        // Since dof_handler is const, we must leave it unchanged.
        const_cast<::Triangulation<DoFHandlerType::dimension,DoFHandlerType::space_dimension> &>(dof_handler.get_triangulation()).load_user_flags (user_flags);
      }

    // Finally, we copy map of sets to
    // map of vectors using the std::vector::assign() function
    std::map< types::global_dof_index, std::vector<typename DoFHandlerType::active_cell_iterator> > dof_to_cell_patches;

    typename std::map<types::global_dof_index, std::set< typename DoFHandlerType::active_cell_iterator> >::iterator
    it = dof_to_set_of_cells_map.begin(),
    it_end = dof_to_set_of_cells_map.end();
    for ( ; it!=it_end; ++it)
      dof_to_cell_patches[it->first].assign( it->second.begin(), it->second.end() );

    return dof_to_cell_patches;
  }

  /*
   * Internally used in collect_periodic_faces
   */
  template <typename CellIterator>
  void
  match_periodic_face_pairs
  (std::set<std::pair<CellIterator, unsigned int> > &pairs1,
   std::set<std::pair<typename identity<CellIterator>::type, unsigned int> > &pairs2,
   const int                                        direction,
   std::vector<PeriodicFacePair<CellIterator> >     &matched_pairs,
   const ::Tensor<1,CellIterator::AccessorType::space_dimension> &offset,
   const FullMatrix<double>                         &matrix)
  {
    static const int space_dim = CellIterator::AccessorType::space_dimension;
    (void)space_dim;
    Assert (0<=direction && direction<space_dim,
            ExcIndexRange (direction, 0, space_dim));

    Assert (pairs1.size() == pairs2.size(),
            ExcMessage ("Unmatched faces on periodic boundaries"));

    unsigned int n_matches = 0;

    // Match with a complexity of O(n^2). This could be improved...
    std::bitset<3> orientation;
    typedef typename std::set
    <std::pair<CellIterator, unsigned int> >::const_iterator PairIterator;
    for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
      {
        for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
          {
            const CellIterator cell1 = it1->first;
            const CellIterator cell2 = it2->first;
            const unsigned int face_idx1 = it1->second;
            const unsigned int face_idx2 = it2->second;
            if (GridTools::orthogonal_equality(orientation,
                                               cell1->face(face_idx1),
                                               cell2->face(face_idx2),
                                               direction, offset,
                                               matrix))
              {
                // We have a match, so insert the matching pairs and
                // remove the matched cell in pairs2 to speed up the
                // matching:
                const PeriodicFacePair<CellIterator> matched_face =
                {
                  {cell1, cell2},
                  {face_idx1, face_idx2},
                  orientation,
                  matrix
                };
                matched_pairs.push_back(matched_face);
                pairs2.erase(it2);
                ++n_matches;
                break;
              }
          }
      }

    //Assure that all faces are matched
    AssertThrow (n_matches == pairs1.size() && pairs2.size() == 0,
                 ExcMessage ("Unmatched faces on periodic boundaries"));
  }



  template <typename MeshType>
  void
  collect_periodic_faces
  (const MeshType                        &mesh,
   const types::boundary_id               b_id,
   const int                              direction,
   std::vector<PeriodicFacePair<typename MeshType::cell_iterator> > &matched_pairs,
   const Tensor<1,MeshType::space_dimension> &offset,
   const FullMatrix<double>              &matrix)
  {
    static const int dim = MeshType::dimension;
    static const int space_dim = MeshType::space_dimension;
    (void)dim;
    (void)space_dim;
    Assert (0<=direction && direction<space_dim,
            ExcIndexRange (direction, 0, space_dim));

    Assert(dim == space_dim,
           ExcNotImplemented());

    // Loop over all cells on the highest level and collect all boundary
    // faces 2*direction and 2*direction*1:

    std::set<std::pair<typename MeshType::cell_iterator, unsigned int> > pairs1;
    std::set<std::pair<typename MeshType::cell_iterator, unsigned int> > pairs2;

    for (typename MeshType::cell_iterator cell = mesh.begin(0);
         cell != mesh.end(0); ++cell)
      {
        const typename MeshType::face_iterator face_1 = cell->face(2*direction);
        const typename MeshType::face_iterator face_2 = cell->face(2*direction+1);

        if (face_1->at_boundary() && face_1->boundary_id() == b_id)
          {
            const std::pair<typename MeshType::cell_iterator, unsigned int> pair1
              = std::make_pair(cell, 2*direction);
            pairs1.insert(pair1);
          }

        if (face_2->at_boundary() && face_2->boundary_id() == b_id)
          {
            const std::pair<typename MeshType::cell_iterator, unsigned int> pair2
              = std::make_pair(cell, 2*direction+1);
            pairs2.insert(pair2);
          }
      }

    Assert (pairs1.size() == pairs2.size(),
            ExcMessage ("Unmatched faces on periodic boundaries"));

    Assert (pairs1.size() > 0,
            ExcMessage("No new periodic face pairs have been found. "
                       "Are you sure that you've selected the correct boundary "
                       "id's and that the coarsest level mesh is colorized?"));

#ifdef DEBUG
    const unsigned int size_old = matched_pairs.size();
#endif

    // and call match_periodic_face_pairs that does the actual matching:
    match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs, offset,
                              matrix);

#ifdef DEBUG
    //check for standard orientation
    const unsigned int size_new = matched_pairs.size();
    for (unsigned int i = size_old; i < size_new; ++i)
      {
        Assert(matched_pairs[i].orientation == 1,
               ExcMessage("Found a face match with non standard orientation. "
                          "This function is only suitable for meshes with cells "
                          "in default orientation"));
      }
#endif
  }



  template <typename MeshType>
  void
  collect_periodic_faces
  (const MeshType                        &mesh,
   const types::boundary_id               b_id1,
   const types::boundary_id               b_id2,
   const int                              direction,
   std::vector<PeriodicFacePair<typename MeshType::cell_iterator> > &matched_pairs,
   const Tensor<1,MeshType::space_dimension> &offset,
   const FullMatrix<double>              &matrix)
  {
    static const int dim = MeshType::dimension;
    static const int space_dim = MeshType::space_dimension;
    (void)dim;
    (void)space_dim;
    Assert (0<=direction && direction<space_dim,
            ExcIndexRange (direction, 0, space_dim));

    // Loop over all cells on the highest level and collect all boundary
    // faces belonging to b_id1 and b_id2:

    std::set<std::pair<typename MeshType::cell_iterator, unsigned int> > pairs1;
    std::set<std::pair<typename MeshType::cell_iterator, unsigned int> > pairs2;

    for (typename MeshType::cell_iterator cell = mesh.begin(0);
         cell != mesh.end(0); ++cell)
      {
        for (unsigned int i = 0; i < GeometryInfo<dim>::faces_per_cell; ++i)
          {
            const typename MeshType::face_iterator face = cell->face(i);
            if (face->at_boundary() && face->boundary_id() == b_id1)
              {
                const std::pair<typename MeshType::cell_iterator, unsigned int> pair1
                  = std::make_pair(cell, i);
                pairs1.insert(pair1);
              }

            if (face->at_boundary() && face->boundary_id() == b_id2)
              {
                const std::pair<typename MeshType::cell_iterator, unsigned int> pair2
                  = std::make_pair(cell, i);
                pairs2.insert(pair2);
              }
          }
      }

    Assert (pairs1.size() == pairs2.size(),
            ExcMessage ("Unmatched faces on periodic boundaries"));

    Assert (pairs1.size() > 0,
            ExcMessage("No new periodic face pairs have been found. "
                       "Are you sure that you've selected the correct boundary "
                       "id's and that the coarsest level mesh is colorized?"));

    // and call match_periodic_face_pairs that does the actual matching:
    match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs, offset,
                              matrix);
  }



  /*
   * Internally used in orthogonal_equality
   *
   * An orthogonal equality test for points:
   *
   * point1 and point2 are considered equal, if
   *   matrix.point1 + offset - point2
   * is parallel to the unit vector in <direction>
   */
  template <int spacedim>
  inline bool orthogonal_equality (const Point<spacedim>    &point1,
                                   const Point<spacedim>    &point2,
                                   const int                 direction,
                                   const Tensor<1,spacedim> &offset,
                                   const FullMatrix<double> &matrix)
  {
    Assert (0<=direction && direction<spacedim,
            ExcIndexRange (direction, 0, spacedim));

    Assert(matrix.m() == matrix.n(), ExcInternalError());

    Point<spacedim> distance;

    if (matrix.m() == spacedim)
      for (int i = 0; i < spacedim; ++i)
        for (int j = 0; j < spacedim; ++j)
          distance(i) += matrix(i,j) * point1(j);
    else
      distance = point1;

    distance += offset - point2;

    for (int i = 0; i < spacedim; ++i)
      {
        // Only compare coordinate-components != direction:
        if (i == direction)
          continue;

        if (fabs(distance(i)) > 1.e-10)
          return false;
      }

    return true;
  }


  /*
   * Internally used in orthogonal_equality
   *
   * A lookup table to transform vertex matchings to orientation flags of
   * the form (face_orientation, face_flip, face_rotation)
   *
   * See the comment on the next function as well as the detailed
   * documentation of make_periodicity_constraints and
   * collect_periodic_faces for details
   */
  template <int dim> struct OrientationLookupTable {};

  template <> struct OrientationLookupTable<1>
  {
    typedef std::array<unsigned int, GeometryInfo<1>::vertices_per_face> MATCH_T;
    static inline std::bitset<3> lookup (const MATCH_T &)
    {
      // The 1D case is trivial
      return 1; // [true ,false,false]
    }
  };

  template <> struct OrientationLookupTable<2>
  {
    typedef std::array<unsigned int, GeometryInfo<2>::vertices_per_face> MATCH_T;
    static inline std::bitset<3> lookup (const MATCH_T &matching)
    {
      // In 2D matching faces (=lines) results in two cases: Either
      // they are aligned or flipped. We store this "line_flip"
      // property somewhat sloppy as "face_flip"
      // (always: face_orientation = true, face_rotation = false)

      static const MATCH_T m_tff = {{ 0, 1 }};
      if (matching == m_tff) return 1;           // [true ,false,false]
      static const MATCH_T m_ttf = {{ 1, 0 }};
      if (matching == m_ttf) return 3;           // [true ,true ,false]
      Assert(false, ExcInternalError());
      // what follows is dead code, but it avoids warnings about the lack
      // of a return value
      return 0;
    }
  };

  template <> struct OrientationLookupTable<3>
  {
    typedef std::array<unsigned int, GeometryInfo<3>::vertices_per_face> MATCH_T;
    static inline std::bitset<3> lookup (const MATCH_T &matching)
    {
      // The full fledged 3D case. *Yay*
      // See the documentation in include/deal.II/base/geometry_info.h
      // as well as the actual implementation in source/grid/tria.cc
      // for more details...

      static const MATCH_T m_tff = {{ 0, 1, 2, 3 }};
      if (matching == m_tff) return 1;                   // [true ,false,false]
      static const MATCH_T m_tft = {{ 1, 3, 0, 2 }};
      if (matching == m_tft) return 5;                   // [true ,false,true ]
      static const MATCH_T m_ttf = {{ 3, 2, 1, 0 }};
      if (matching == m_ttf) return 3;                   // [true ,true ,false]
      static const MATCH_T m_ttt = {{ 2, 0, 3, 1 }};
      if (matching == m_ttt) return 7;                   // [true ,true ,true ]
      static const MATCH_T m_fff = {{ 0, 2, 1, 3 }};
      if (matching == m_fff) return 0;                   // [false,false,false]
      static const MATCH_T m_fft = {{ 2, 3, 0, 1 }};
      if (matching == m_fft) return 4;                   // [false,false,true ]
      static const MATCH_T m_ftf = {{ 3, 1, 2, 0 }};
      if (matching == m_ftf) return 2;                   // [false,true ,false]
      static const MATCH_T m_ftt = {{ 1, 0, 3, 2 }};
      if (matching == m_ftt) return 6;                   // [false,true ,true ]
      Assert(false, ExcInternalError());
      // what follows is dead code, but it avoids warnings about the lack
      // of a return value
      return 0;
    }
  };



  template <typename FaceIterator>
  inline bool
  orthogonal_equality (std::bitset<3>     &orientation,
                       const FaceIterator &face1,
                       const FaceIterator &face2,
                       const int          direction,
                       const Tensor<1,FaceIterator::AccessorType::space_dimension> &offset,
                       const FullMatrix<double> &matrix)
  {
    Assert(matrix.m() == matrix.n(),
           ExcMessage("The supplied matrix must be a square matrix"));

    static const int dim = FaceIterator::AccessorType::dimension;

    // Do a full matching of the face vertices:

    std::
    array<unsigned int, GeometryInfo<dim>::vertices_per_face> matching;

    std::set<unsigned int> face2_vertices;
    for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face; ++i)
      face2_vertices.insert(i);

    for (unsigned int i = 0; i < GeometryInfo<dim>::vertices_per_face; ++i)
      {
        for (std::set<unsigned int>::iterator it = face2_vertices.begin();
             it != face2_vertices.end();
             ++it)
          {
            if (orthogonal_equality(face1->vertex(i),face2->vertex(*it),
                                    direction, offset, matrix))
              {
                matching[i] = *it;
                face2_vertices.erase(it);
                break; // jump out of the innermost loop
              }
          }
      }

    // And finally, a lookup to determine the ordering bitmask:
    if (face2_vertices.empty())
      orientation = OrientationLookupTable<dim>::lookup(matching);

    return face2_vertices.empty();
  }



  template <typename FaceIterator>
  inline bool
  orthogonal_equality (const FaceIterator &face1,
                       const FaceIterator &face2,
                       const int          direction,
                       const Tensor<1,FaceIterator::AccessorType::space_dimension> &offset,
                       const FullMatrix<double> &matrix)
  {
    // Call the function above with a dummy orientation array
    std::bitset<3> dummy;
    return orthogonal_equality (dummy, face1, face2, direction, offset, matrix);
  }



}


#include "grid_tools_dof_handlers.inst"


DEAL_II_NAMESPACE_CLOSE