From: Spencer Patty <srobertp@gmail.com>
Date: Thu, 4 Feb 2016 17:40:39 +0000 (-0600)
Subject: Merge remote-tracking branch 'origin/master' into dof_support_patches
X-Git-Tag: v8.4.0-rc2~19^2~8
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c5966d2c0dfd2ef0b7ef05ca910589ecff521411;p=dealii.git

Merge remote-tracking branch 'origin/master' into dof_support_patches
---

c5966d2c0dfd2ef0b7ef05ca910589ecff521411
diff --cc include/deal.II/grid/grid_tools.h
index 1daf92ad7a,13bca31ef3..79920d8459
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@@ -974,41 -978,107 +978,142 @@@ namespace GridTool
    get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
  
  
+   /**
+    * This function takes a vector of active cells (hereafter named @p
+    * patch_cells) as input argument, and returns a vector of their
+    * parent cells with the coarsest common level of refinement. In
+    * other words, find that set of cells living at the same refinement
+    * level so that all cells in the input vector are children of the
+    * cells in the set, or are in the set itself.
+    *
+    * @tparam Container In C++, the compiler can not determine the type
+    * of <code>Container</code> from the function call. You need to
+    * specify it as an explicit template argument following the
+    * function name. This type has to satisfy the requirements of a
+    * mesh container (see @ref GlossMeshAsAContainer).
+    *
+    * @param[in] patch_cells A vector of active cells for which
+    *   this function finds the parents at the coarsest common
+    *   level. This vector of cells typically results from
+    *   calling the function GridTools::get_patch_around_cell().
+    * @return A list of cells with the coarsest common level of
+    *   refinement of the input cells.
+    *
+    * @author Arezou Ghesmati, Wolfgang Bangerth, 2015
+    */
+   template <class Container>
+   std::vector<typename Container::cell_iterator>
+   get_cells_at_coarsest_common_level(const std::vector<typename Container::active_cell_iterator> &patch_cells);
+ 
+   /**
+    * This function constructs a Triangulation (named @p
+    * local_triangulation) from a given vector of active cells. This
+    * vector (which we think of the cells corresponding to a "patch")
+    * contains active cells that are part of an existing global
+    * Triangulation. The goal of this function is to build a local
+    * Triangulation that contains only the active cells given in
+    * @p patch (and potentially a minimum number of additional cells
+    * required to form a valid Triangulation).
+    * The function also returns a map that allows to identify the cells in
+    * the output Triangulation and corresponding cells in the input
+    * list.
+    *
+    * The operation implemented by this function is frequently used in
+    * the definition of error estimators that need to solve "local"
+    * problems on each cell and its neighbors. A similar construction is
+    * necessary in the definition of the Clement interpolation operator
+    * in which one needs to solve a local problem on all cells within
+    * the support of a shape function. This function then builds a
+    * complete Triangulation from a list of cells that make up such a
+    * patch; one can then later attach a DoFHandler to such a
+    * Triangulation.
+    *
+    * If the list of input cells contains only cells at the same
+    * refinement level, then the output Triangulation simply consists
+    * of a Triangulation containing only exactly these patch cells. On
+    * the other hand, if the input cells live on different refinement
+    * levels, i.e., the Triangulation of which they are part is
+    * adaptively refined, then the construction of the output
+    * Triangulation is not so simple because the coarsest level of a
+    * Triangulation can not contain hanging nodes. Rather, we first
+    * have to find the common refinement level of all input cells,
+    * along with their common parents (see
+    * GridTools::get_cells_at_coarsest_common_level()), build a
+    * Triangulation from those, and then adaptively refine it so that
+    * the input cells all also exist in the output Triangulation.
+    *
+    * A consequence of this procedure is that that output Triangulation
+    * may contain more active cells than the ones that exist in the
+    * input vector.  On the other hand, one typically wants to solve
+    * the local problem not on the entire output Triangulation, but
+    * only on those cells of it that correspond to cells in the input
+    * list.  In this case, a user typically wants to assign degrees of
+    * freedom only on cells that are part of the "patch", and somehow
+    * ignore those excessive cells. The current function supports this
+    * common requirement by setting the user flag for the cells in the
+    * output Triangulation that match with cells in the input
+    * list. Cells which are not part of the original patch will not
+    * have their @p user_flag set; we can then avoid assigning degrees of
+    * freedom using the FE_Nothing<dim> element.
+    *
+    * @tparam Container In C++, the compiler can not determine the type
+    * of <code>Container</code> from the function call. You need to
+    * specify it as an explicit template argument following the
+    * function name. This type that satisfies the requirements of a
+    * mesh container (see @ref GlossMeshAsAContainer).
+    *
+    * @param[in] patch A vector of active cells from a common triangulation.
+    *  These cells may or may not all be at the same refinement level.
+    * @param[out] local_triangulation A triangulation whose active cells
+    *  correspond to the given vector of active cells in @p patch.
+    * @param[out] patch_to_global_tria_map A map between the local triangulation
+    * which is built as explained above, and the cell iterators in the input list.
+    *
+    *  @author Arezou Ghesmati, Wolfgang Bangerth, 2015
+    */
+   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);
+ 
 +  /**
 +   * This function runs through the degrees of freedom defined by the
 +   * Container and for each dof constructs a vector of active_cell_iterators
 +   * representing the cells of support of the nodal basis element
 +   * at that degree of freedom. This function was designed for the
 +   * implementation of local projections, for instance the Clement interpolant,
 +   * in conjunction with other local patch functions like
 +   * GridTools::build_triangulation_from_patch.
 +   *
 +   * Containers built on top of Triangulation<dim> or
 +   * parallel:distributed::Triangulation<dim> are supported and handled
 +   * appropriately.
 +   *
 +   * It is necessary that the finite element underlying the Container used has
 +   * degrees of freedom on faces (2d or 3d) and lines (in 3d).  This unfortunately
 +   * precludes the FE_DGQ<dim> finite element.  Likewise, the finite element
 +   * must have nodal basis elements for this implementation to make sense.
 +   *
 +   * @tparam Container The Container could be a DoFHandler<dim> or hp::DoFHandler<dim>.
 +   * @param[in] dof_handler The Container which could be built on a Triangulation<dim>
 +   * or a parallel:distributed::Triangulation<dim> and should be using a nodal
 +   * finite element with degrees of freedom defined on faces (2d or 3d) and
 +   * lines (3d).
 +   * @param[out] dof_to_support_patch_map A map from the global_dof_index of dofs on locally relevant cells
 +   * to vectors containing Container::active_cell_iterators of
 +   * cells in support of basis function at that degree of freedom.
 +   *
 +   *  @author Spencer Patty, 2016
 +   *
 +   */
 +  template <class Container>
 +  std::map< types::global_dof_index,std::vector<typename Container::active_cell_iterator> >
 +  get_dof_to_support_patch_map(Container& dof_handler);
 +
 +
    /*@}*/
    /**
     * @name Lower-dimensional meshes for parts of higher-dimensional meshes
diff --cc source/grid/grid_tools.cc
index b003cea190,48bdeb3c6a..d06c53e1dc
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@@ -2872,256 -2872,182 +2872,431 @@@ next_cell
  
  
  
+   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)
+       {
+         bool repeat_vertex;
+         for (unsigned int j=0;  j< GeometryInfo<Container::dimension>::vertices_per_cell; ++j)
+           {
+             Point<Container::space_dimension> position=(*uniform_cell)->vertex (j);
+             repeat_vertex=false;
+ 
+             for (unsigned int m=0; m<i; ++m)
+               {
+                 if (position == vertices[m])
+                   {
+                     repeat_vertex=true;
+                     cells[k].vertices[j]=m;
+                     break;
+                   }
+               }
+             if (repeat_vertex==false)
+               {
+                 vertices.push_back(position);
+                 cells[k].vertices[j]=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)
+                 {
+                   if (patch_to_global_tria_map_tmp[active_tria_cell]==patch[i])
+                     {
+                       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)));
+ 
+                                 Assert(cell->child(c)->center().distance( patch_to_global_tria_map_tmp[cell]->child(c)->center())
+                                        <=1e-15*cell->child(c)->diameter(),
+                                        ExcInternalError());
+                               }
+                           }
+                         // 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);
+     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 Container>
 +  std::map< types::global_dof_index,std::vector<typename Container::active_cell_iterator> >
 +  get_dof_to_support_patch_map(Container& dof_handler)
 +  {
 +
 +    //  TODO: Add Assert( fe is not dg)
 +
 +    // 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 Container::active_cell_iterator> > dof_to_set_of_cells_map;
 +
 +    std::vector<types::global_dof_index> local_dof_index;
 +    std::vector<types::global_dof_index> local_face_dof_index;
 +    std::vector<types::global_dof_index> local_line_dof_index;
 +
 +    // in 3d, we need pointers from active lines to the
 +    // active parent lines, so we construct it as needed.
 +    std::map<typename Container::active_line_iterator, typename Container::line_iterator > lines_to_parent_lines_map;
 +    if (Container::dimension == 3)
 +    {
 +      typename Container::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<Container::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 Container::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_index.resize(n_dofs_per_cell);
 +
 +        // Take care of adding cell pointer to each
 +        // dofs that exists on cell.
 +        cell->get_dof_indices(local_dof_index);
 +        for (unsigned int i=0; i< n_dofs_per_cell; ++i )
 +          dof_to_set_of_cells_map[local_dof_index[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<Container::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_index.resize(n_dofs_per_face);
 +
 +              cell->face(f)->child(c)->get_dof_indices(local_face_dof_index);
 +              for (unsigned int i=0; i< n_dofs_per_face; ++i )
 +                dof_to_set_of_cells_map[local_face_dof_index[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_index.resize(n_dofs_per_face);
 +
 +            cell->neighbor(f)->face(face_no)->get_dof_indices(local_face_dof_index);
 +            for (unsigned int i=0; i< n_dofs_per_face; ++i )
 +              dof_to_set_of_cells_map[local_face_dof_index[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_index.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_index);
 +                for (unsigned int i=0; i<n_dofs_per_face; ++i )
 +                  dof_to_set_of_cells_map[local_face_dof_index[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 (Container::dimension == 3)
 +        {
 +          for (unsigned int l=0; l<GeometryInfo<Container::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_index.resize(n_dofs_per_line);
 +
 +                cell->line(l)->child(c)->get_dof_indices(local_line_dof_index);
 +                for (unsigned int i=0; i<n_dofs_per_line; ++i )
 +                  dof_to_set_of_cells_map[local_line_dof_index[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 Container::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_index.resize(n_dofs_per_line);
 +
 +              parent_line->get_dof_indices(local_line_dof_index);
 +              for (unsigned int i=0; i<n_dofs_per_line; ++i )
 +                dof_to_set_of_cells_map[local_line_dof_index[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_index.resize(n_dofs_per_line);
 +
 +                parent_line->child(c)->get_dof_indices(local_line_dof_index);
 +                for (unsigned int i=0; i<n_dofs_per_line; ++i )
 +                  dof_to_set_of_cells_map[local_line_dof_index[i]].insert(cell);
 +              }
 +
 +              // clear up user flags set from earlier
 +              // denoting that a 3d line has an active
 +              // parent and so an entry exists in the
 +              // lines_to_parent_lines_map map for that
 +              // line pointing to it's parent.
 +              cell->line(l)->clear_user_flag();
 +            }
 +          } // for lines l
 +        }// if Container::dimension == 3
 +      }// if cell->is_locally_owned()
 +    }// for cells
 +
 +
 +    // Finally, we copy map of sets to
 +    // map of vectors using assign()
 +    std::map< types::global_dof_index, std::vector<typename Container::active_cell_iterator> > dof_to_cell_patches;
 +
 +    typename std::map<types::global_dof_index, std::set< typename Container::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
     *