]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Split the grid tools instantiations.
authorDavid Wells <wellsd2@rpi.edu>
Sun, 28 Jan 2018 23:49:55 +0000 (18:49 -0500)
committerDavid Wells <wellsd2@rpi.edu>
Sat, 3 Feb 2018 18:22:42 +0000 (13:22 -0500)
This commit splits the grid tools instantiation files in two: one file deals
with instantiations for DoFHandlers and Triangulations and the other deals with
instantiations of just Triangulations. Here is the timing information:

before this patch:
debug:
grid tools: 74s, 3.6 GB

release:
grid tools: 97s 3.9 GB

after this patch:
debug:
grid tools: 31s, 1.8 GB
grid tools dof handlers: 48s, 2.4 GB

release:
grid tools: 46s, 1.9 GB
grid tools dof handlers: 65s, 2.4 GB

This lowers the high water mark for memory usage with only a small impact on the
total compilation time.

source/grid/CMakeLists.txt
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in
source/grid/grid_tools_dof_handlers.cc [new file with mode: 0644]
source/grid/grid_tools_dof_handlers.inst.in [new file with mode: 0644]

index a5888e4730a34c556975d352a477bddf5a1b7c7c..a1e4c1e165d0023b8b6e63d71d78119dd42f3906 100644 (file)
@@ -37,6 +37,7 @@ SET(_unity_include_src
 SET(_separate_src
   grid_reordering.cc
   grid_tools.cc
+  grid_tools_dof_handlers.cc
   tria.cc
   )
 
@@ -56,6 +57,7 @@ SET(_inst
   grid_out.inst.in
   grid_refinement.inst.in
   grid_tools.inst.in
+  grid_tools_dof_handlers.inst.in
   grid_tools_cache.inst.in
   intergrid_map.inst.in
   manifold.inst.in
index 43e7f53cbff6bec2824e0405e217226a474270f8..3dd3083319331ba94945c36983fbe9299d597fc1 100644 (file)
@@ -44,7 +44,6 @@
 #include <deal.II/fe/mapping_q1.h>
 #include <deal.II/fe/mapping_q.h>
 #include <deal.II/fe/fe_values.h>
-#include <deal.II/hp/mapping_collection.h>
 #include <deal.II/numerics/matrix_tools.h>
 #include <deal.II/lac/sparsity_tools.h>
 
@@ -1068,477 +1067,6 @@ 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 dealii::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 dealii::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type> adjacent_cells;
-
-    typename dealii::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 dealii::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 dealii::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type> &searched_cells,
-     std::set<typename dealii::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 dealii::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 dealii::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 dealii::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 dealii::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 <int dim,int spacedim>
   std::vector<std::vector<Tensor<1,spacedim> > >
   vertex_to_cell_centers_directions(const Triangulation<dim,spacedim> &mesh,
@@ -1735,465 +1263,6 @@ next_cell:
 
 
 
-  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
-                  }
-              }
-            //udpate 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;
-  }
-
-
-  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_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;
-  }
-
-
-
-  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 );
-  }
-
-
-
   namespace internal
   {
     namespace BoundingBoxPredicate
@@ -3150,143 +2219,9 @@ next_cell:
 
 
 
-  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>
-  double
-  minimal_cell_diameter (const Triangulation<dim, spacedim> &triangulation)
+  template <int dim, int spacedim>
+  double
+  minimal_cell_diameter (const Triangulation<dim, spacedim> &triangulation)
   {
     double min_diameter = triangulation.begin_active()->diameter();
     for (typename Triangulation<dim, spacedim>::active_cell_iterator
@@ -3790,920 +2725,6 @@ next_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 prefered 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<dealii::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<dealii::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 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);
-  }
-
-
-
-  /*
-   * 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 dealii::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_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);
-  }
-
-
-
-  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 <int dim, int spacedim>
   void copy_boundary_to_manifold_id(Triangulation<dim, spacedim> &tria,
                                     const bool reset_boundary_ids)
index 51d4df9a32c8027d4cde74ef6b3efecb51c68004..22d8a2afdbd55be23cba8ecbc3cb208e586103aa 100644 (file)
@@ -39,7 +39,8 @@ for (X : TRIANGULATIONS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimensio
 }
 
 
-for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
+// now also instantiate a few additional functions for parallel::distributed::Triangulation
+for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
 {
 
 #if deal_II_dimension <= deal_II_space_dimension
@@ -47,122 +48,38 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II
 
     template
     unsigned int
-    find_closest_vertex (const X &,
-                         const Point<deal_II_space_dimension> &,
-                         const std::vector<bool> &);
-
-    template
-    unsigned int
-    find_closest_vertex (const Mapping<deal_II_dimension, deal_II_space_dimension>&,
-                         const X &,
-                         const Point<deal_II_space_dimension> &,
-                         const std::vector<bool> &);
-
-    template
-    std::vector<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type>
-    find_cells_adjacent_to_vertex(const X &, const unsigned int);
-
-    template
-    dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type
-    find_active_cell_around_point (const X &, const Point<deal_II_space_dimension> &, const std::vector<bool> &);
-
-    template
-    std::pair<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type, Point<deal_II_dimension> >
-    find_active_cell_around_point (const Mapping<deal_II_dimension, deal_II_space_dimension> &,
-                                   const X &,
-                                   const Point<deal_II_space_dimension> &,
-                                   const std::vector<bool> &);
-
-    template
-    std::vector<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type>
-    compute_active_cell_halo_layer (const X &,
-                                    const std::function<bool (const dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type&)> &);
-
-    template
-    std::vector<X::cell_iterator>
-    compute_cell_halo_layer_on_level (const X &,
-                                      const std::function<bool (const X::cell_iterator&)> &,
-                                      const unsigned int);
-
-    template
-    std::vector<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type>
-    compute_ghost_cell_halo_layer (const X &);
-
-
-    template
-    std::vector<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type>
-    compute_active_cell_layer_within_distance (const X &,
-            const std::function<bool (const dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type&)> &,
-            const double);
-
+    find_closest_vertex_of_cell<deal_II_dimension,deal_II_space_dimension>
+    (const typename Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator &,
+     const Point<deal_II_space_dimension> &);
 
     template
-    std::vector<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type>
-    compute_ghost_cell_layer_within_distance (const X &, const double);
-
+    std::map<unsigned int,types::global_vertex_index>
+    compute_local_to_global_vertex_index_map(const parallel::distributed::Triangulation<deal_II_dimension,deal_II_space_dimension> &triangulation);
 
     template
-    std::pair< Point<X::space_dimension>, Point<X::space_dimension> >
-    compute_bounding_box (const X &,
-                          const std::function<bool (const dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type&)> &);
-
+    std::map<unsigned int, Point<deal_II_space_dimension> >
+    extract_used_vertices(const Triangulation<deal_II_dimension,deal_II_space_dimension>&mesh,
+                          const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping);
 
     template
-    std::list<std::pair<X::cell_iterator, X::cell_iterator> >
-    get_finest_common_cells (const X &mesh_1,
-                             const X &mesh_2);
-
+    std::pair<typename Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator,
+              Point<deal_II_dimension> >
+    find_active_cell_around_point(const Cache<deal_II_dimension,deal_II_space_dimension>&,
+                                  const Point<deal_II_space_dimension> &,
+                                  const typename Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator &,
+                                  const std::vector<bool>  &);
 
     template
-    bool
-    have_same_coarse_mesh (const X &mesh_1,
-                           const X &mesh_2);
+    std::tuple< std::vector< typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator >,
+                std::vector< std::vector< Point< deal_II_dimension > > >, std::vector< std::vector< unsigned int > > >
+    compute_point_locations(const Cache< deal_II_dimension, deal_II_space_dimension > &,
+                            const std::vector< Point< deal_II_space_dimension > > &,
+                            const typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator &);
     \}
 
 #endif
 }
 
-// now also instantiate a few additional functions for parallel::distributed::Triangulation
-for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
-{
-
-#if deal_II_dimension <= deal_II_space_dimension
-    namespace GridTools \{
-
-        template
-        unsigned int
-        find_closest_vertex_of_cell<deal_II_dimension,deal_II_space_dimension>
-        (const typename Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator &,
-         const Point<deal_II_space_dimension> &);
-
-        template
-        std::map<unsigned int,types::global_vertex_index>
-        compute_local_to_global_vertex_index_map(const parallel::distributed::Triangulation<deal_II_dimension,deal_II_space_dimension> &triangulation);
-
-        template
-        std::map<unsigned int, Point<deal_II_space_dimension> >
-        extract_used_vertices(const Triangulation<deal_II_dimension,deal_II_space_dimension>&mesh,
-                              const Mapping<deal_II_dimension,deal_II_space_dimension> &mapping);
-
-        template
-        std::pair<typename Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator,
-                  Point<deal_II_dimension> >
-        find_active_cell_around_point(const Cache<deal_II_dimension,deal_II_space_dimension>&,
-                                      const Point<deal_II_space_dimension> &,
-                                      const typename Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator &,
-                                      const std::vector<bool>  &);
-
-        template
-        std::tuple< std::vector< typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator >,
-                    std::vector< std::vector< Point< deal_II_dimension > > >, std::vector< std::vector< unsigned int > > >
-        compute_point_locations(const Cache< deal_II_dimension, deal_II_space_dimension > &,
-                                const std::vector< Point< deal_II_space_dimension > > &,
-                                const typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator &);
-                                                                                                      \}
-
-#endif
-}
-
 
 
 for (deal_II_space_dimension : SPACE_DIMENSIONS)
@@ -284,14 +201,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
     template
     void partition_multigrid_levels (Triangulation<deal_II_dimension, deal_II_space_dimension> &);
 
-    template
-    std::pair<hp::DoFHandler<deal_II_dimension, deal_II_space_dimension>::active_cell_iterator,
-              Point<deal_II_dimension> >
-    find_active_cell_around_point
-    (const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension> &,
-     const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
-     const Point<deal_II_space_dimension> &);
-
     template
     void get_subdomain_association (const Triangulation<deal_II_dimension, deal_II_space_dimension>  &,
                                     std::vector<types::subdomain_id> &);
@@ -351,141 +260,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
 
 
 
-}
-
-for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS; Container : DOFHANDLER_TEMPLATES)
-{
-#if deal_II_dimension <= deal_II_space_dimension
-    namespace GridTools \{
-
-    template
-    std::map< types::global_dof_index,std::vector<Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator> >
-    get_dof_to_support_patch_map<Container<deal_II_dimension,deal_II_space_dimension> >
-    (Container<deal_II_dimension,deal_II_space_dimension> &dof_handler);
-
-    \}
-#endif
-}
-
-
-for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS; Container : TRIANGULATION_AND_DOFHANDLER_TEMPLATES)
-{
-#if deal_II_dimension <= deal_II_space_dimension
-    namespace GridTools \{
-
-    template
-    std::vector<Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator>
-    get_patch_around_cell<Container<deal_II_dimension,deal_II_space_dimension> >
-    (const Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator &cell);
-
-    template
-    std::vector< Container<deal_II_dimension,deal_II_space_dimension>::cell_iterator>
-    get_cells_at_coarsest_common_level <Container<deal_II_dimension,deal_II_space_dimension> > (
-        const std::vector< Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator> & patch_cells);
-
-    template
-    void build_triangulation_from_patch <Container<deal_II_dimension,deal_II_space_dimension> > (
-        const std::vector<Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator>  &patch,
-        Triangulation<Container<deal_II_dimension,deal_II_space_dimension>::dimension,Container<deal_II_dimension,deal_II_space_dimension>::space_dimension> &local_triangulation,
-        std::map<Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator,
-        Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator > &patch_to_global_tria_map);
-
-    \}
-#endif
-}
-
-
-// instantiate the following functions only for the "sequential" containers. this
-// is a misnomer here, however: the point is simply that we only instantiate
-// these functions for certain *iterator* types, and the iterator types are
-// the same for sequential and parallel containers; consequently, we get duplicate
-// instantiation errors if we instantiate for *all* container types, rather than
-// only the sequential ones
-for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
-{
-#if deal_II_dimension <= deal_II_space_dimension
-    namespace GridTools \{
-
-    template
-    bool orthogonal_equality<X::active_face_iterator> (std::bitset<3> &,
-            const X::active_face_iterator&,
-            const X::active_face_iterator&,
-            const int,
-            const Tensor<1,deal_II_space_dimension> &,
-            const FullMatrix<double> &);
-
-    template
-    bool orthogonal_equality<X::face_iterator> (std::bitset<3> &,
-            const X::face_iterator&,
-            const X::face_iterator&,
-            const int,
-            const Tensor<1,deal_II_space_dimension> &,
-            const FullMatrix<double> &);
-
-    template
-    bool orthogonal_equality<X::active_face_iterator> (const X::active_face_iterator&,
-            const X::active_face_iterator&,
-            const int,
-            const Tensor<1,deal_II_space_dimension> &,
-            const FullMatrix<double> &);
-
-    template
-    bool orthogonal_equality<X::face_iterator> (const X::face_iterator&,
-            const X::face_iterator&,
-            const int,
-            const Tensor<1,deal_II_space_dimension> &,
-            const FullMatrix<double> &);
-
-    template
-    void collect_periodic_faces<X> (const X &,
-                                    const types::boundary_id,
-                                    const types::boundary_id,
-                                    const int,
-                                    std::vector<PeriodicFacePair<X::cell_iterator> > &,
-                                    const Tensor<1,X::space_dimension> &,
-                                    const FullMatrix<double> &);
-
-    template
-    void collect_periodic_faces<X> (const X &,
-                                    const types::boundary_id,
-                                    const int,
-                                    std::vector<PeriodicFacePair<X::cell_iterator> > &,
-                                    const Tensor<1,X::space_dimension> &,
-                                    const FullMatrix<double> &);
-
-    \}
-#endif
-}
-
-for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
-{
-#if deal_II_dimension <= deal_II_space_dimension
-#if deal_II_dimension >= 2
-
-    namespace GridTools \{
-    template
-    void
-    collect_periodic_faces<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> >
-    (const parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> &,
-     const types::boundary_id,
-     const types::boundary_id,
-     const int,
-     std::vector<PeriodicFacePair<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator> > &,
-     const Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &,
-     const FullMatrix<double> &);
-
-    template
-    void
-    collect_periodic_faces<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> >
-    (const parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> &,
-     const types::boundary_id,
-     const int,
-     std::vector<PeriodicFacePair<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator> > &,
-     const Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &,
-     const FullMatrix<double> &);
-    \}
-#endif
-#endif
 }
 
 for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc
new file mode 100644 (file)
index 0000000..2caf255
--- /dev/null
@@ -0,0 +1,2050 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 2017 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/quadrature_lib.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_memory.h>
+#include <deal.II/lac/filtered_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/distributed/shared_tria.h>
+#include <deal.II/distributed/tria_base.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_boundary.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+#include <deal.II/grid/grid_reordering.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/hp/mapping_collection.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/lac/sparsity_tools.h>
+
+#include <boost/random/uniform_real_distribution.hpp>
+#include <boost/random/mersenne_twister.hpp>
+
+#include <array>
+#include <cmath>
+#include <numeric>
+#include <list>
+#include <set>
+#include <tuple>
+
+
+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 dealii::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 dealii::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type> adjacent_cells;
+
+    typename dealii::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 dealii::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 dealii::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type> &searched_cells,
+     std::set<typename dealii::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 dealii::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 dealii::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 dealii::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 dealii::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
+                  }
+              }
+            //udpate 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 prefered 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<dealii::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<dealii::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 dealii::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
diff --git a/source/grid/grid_tools_dof_handlers.inst.in b/source/grid/grid_tools_dof_handlers.inst.in
new file mode 100644 (file)
index 0000000..1444b26
--- /dev/null
@@ -0,0 +1,242 @@
+
+for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
+{
+
+#if deal_II_dimension <= deal_II_space_dimension
+    namespace GridTools \{
+
+    template
+    unsigned int
+    find_closest_vertex (const X &,
+                         const Point<deal_II_space_dimension> &,
+                         const std::vector<bool> &);
+
+    template
+    unsigned int
+    find_closest_vertex (const Mapping<deal_II_dimension, deal_II_space_dimension>&,
+                         const X &,
+                         const Point<deal_II_space_dimension> &,
+                         const std::vector<bool> &);
+
+    template
+    std::vector<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type>
+    find_cells_adjacent_to_vertex(const X &,
+                                  const unsigned int);
+
+    template
+    dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type
+    find_active_cell_around_point (const X &,
+                                   const Point<deal_II_space_dimension> &,
+                                   const std::vector<bool> &);
+
+    template
+    std::pair<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type, Point<deal_II_dimension> >
+    find_active_cell_around_point (const Mapping<deal_II_dimension, deal_II_space_dimension> &,
+                                   const X &,
+                                   const Point<deal_II_space_dimension> &,
+                                   const std::vector<bool> &);
+
+    template
+    std::vector<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type>
+    compute_active_cell_halo_layer (const X &,
+                                    const std::function<bool (const dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type&)> &);
+
+    template
+    std::vector<X::cell_iterator>
+    compute_cell_halo_layer_on_level (const X &,
+                                      const std::function<bool (const X::cell_iterator&)> &,
+                                      const unsigned int);
+
+    template
+    std::vector<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type>
+    compute_ghost_cell_halo_layer (const X &);
+
+
+    template
+    std::vector<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type>
+    compute_active_cell_layer_within_distance (const X &,
+            const std::function<bool (const dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type&)> &,
+            const double);
+
+
+    template
+    std::vector<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type>
+    compute_ghost_cell_layer_within_distance (const X &, const double);
+
+
+    template
+    std::pair< Point<X::space_dimension>, Point<X::space_dimension> >
+    compute_bounding_box (const X &,
+                          const std::function<bool (const dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type&)> &);
+
+
+    template
+    std::list<std::pair<X::cell_iterator, X::cell_iterator> >
+    get_finest_common_cells (const X &mesh_1,
+                             const X &mesh_2);
+
+
+    template
+    bool
+    have_same_coarse_mesh (const X &mesh_1,
+                           const X &mesh_2);
+    \}
+
+#endif
+}
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+{
+#if deal_II_dimension <= deal_II_space_dimension
+    namespace GridTools \{
+
+        template
+        std::pair<hp::DoFHandler<deal_II_dimension, deal_II_space_dimension>::active_cell_iterator,
+                  Point<deal_II_dimension> >
+        find_active_cell_around_point
+        (const hp::MappingCollection<deal_II_dimension, deal_II_space_dimension> &,
+         const hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> &,
+         const Point<deal_II_space_dimension> &);
+
+                       \}
+#endif
+}
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS; Container : TRIANGULATION_AND_DOFHANDLER_TEMPLATES)
+{
+#if deal_II_dimension <= deal_II_space_dimension
+    namespace GridTools \{
+
+    template
+    std::vector<Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator>
+    get_patch_around_cell<Container<deal_II_dimension,deal_II_space_dimension> >
+    (const Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator &cell);
+
+    template
+    std::vector< Container<deal_II_dimension,deal_II_space_dimension>::cell_iterator>
+    get_cells_at_coarsest_common_level <Container<deal_II_dimension,deal_II_space_dimension> > (
+        const std::vector< Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator> & patch_cells);
+
+    template
+    void build_triangulation_from_patch <Container<deal_II_dimension,deal_II_space_dimension> > (
+        const std::vector<Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator>  &patch,
+        Triangulation<Container<deal_II_dimension,deal_II_space_dimension>::dimension,Container<deal_II_dimension,deal_II_space_dimension>::space_dimension> &local_triangulation,
+        std::map<Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator,
+        Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator > &patch_to_global_tria_map);
+
+    \}
+#endif
+}
+
+
+
+for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS; Container : DOFHANDLER_TEMPLATES)
+{
+#if deal_II_dimension <= deal_II_space_dimension
+    namespace GridTools \{
+
+    template
+    std::map< types::global_dof_index,std::vector<Container<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator> >
+    get_dof_to_support_patch_map<Container<deal_II_dimension,deal_II_space_dimension> >
+    (Container<deal_II_dimension,deal_II_space_dimension> &dof_handler);
+
+    \}
+#endif
+}
+
+
+
+// instantiate the following functions only for the "sequential" containers. this
+// is a misnomer here, however: the point is simply that we only instantiate
+// these functions for certain *iterator* types, and the iterator types are
+// the same for sequential and parallel containers; consequently, we get duplicate
+// instantiation errors if we instantiate for *all* container types, rather than
+// only the sequential ones
+for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
+{
+#if deal_II_dimension <= deal_II_space_dimension
+    namespace GridTools \{
+
+    template
+    bool orthogonal_equality<X::active_face_iterator> (std::bitset<3> &,
+            const X::active_face_iterator&,
+            const X::active_face_iterator&,
+            const int,
+            const Tensor<1,deal_II_space_dimension> &,
+            const FullMatrix<double> &);
+
+    template
+    bool orthogonal_equality<X::face_iterator> (std::bitset<3> &,
+            const X::face_iterator&,
+            const X::face_iterator&,
+            const int,
+            const Tensor<1,deal_II_space_dimension> &,
+            const FullMatrix<double> &);
+
+    template
+    bool orthogonal_equality<X::active_face_iterator> (const X::active_face_iterator&,
+            const X::active_face_iterator&,
+            const int,
+            const Tensor<1,deal_II_space_dimension> &,
+            const FullMatrix<double> &);
+
+    template
+    bool orthogonal_equality<X::face_iterator> (const X::face_iterator&,
+            const X::face_iterator&,
+            const int,
+            const Tensor<1,deal_II_space_dimension> &,
+            const FullMatrix<double> &);
+
+    template
+    void collect_periodic_faces<X> (const X &,
+                                    const types::boundary_id,
+                                    const types::boundary_id,
+                                    const int,
+                                    std::vector<PeriodicFacePair<X::cell_iterator> > &,
+                                    const Tensor<1,X::space_dimension> &,
+                                    const FullMatrix<double> &);
+
+    template
+    void collect_periodic_faces<X> (const X &,
+                                    const types::boundary_id,
+                                    const int,
+                                    std::vector<PeriodicFacePair<X::cell_iterator> > &,
+                                    const Tensor<1,X::space_dimension> &,
+                                    const FullMatrix<double> &);
+
+    \}
+#endif
+}
+
+// TODO the text above the last instantiation block implies that this should not be necessary... is it?
+for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
+{
+#if deal_II_dimension <= deal_II_space_dimension
+#if deal_II_dimension >= 2
+
+    namespace GridTools \{
+    template
+    void
+    collect_periodic_faces<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> >
+    (const parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> &,
+     const types::boundary_id,
+     const types::boundary_id,
+     const int,
+     std::vector<PeriodicFacePair<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator> > &,
+     const Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &,
+     const FullMatrix<double> &);
+
+    template
+    void
+    collect_periodic_faces<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> >
+    (const parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension> &,
+     const types::boundary_id,
+     const int,
+     std::vector<PeriodicFacePair<parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::cell_iterator> > &,
+     const Tensor<1,parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>::space_dimension> &,
+     const FullMatrix<double> &);
+    \}
+#endif
+#endif
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.