]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add the get_dof_to_support_patch_map() function to grid_tools and instantiations.
authorSpencer Patty <srobertp@gmail.com>
Wed, 3 Feb 2016 19:21:50 +0000 (13:21 -0600)
committerSpencer Patty <srobertp@gmail.com>
Wed, 3 Feb 2016 19:21:50 +0000 (13:21 -0600)
include/deal.II/grid/grid_tools.h
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in

index 268c96f378b686b137669ec21c19c223b9ec3c13..1daf92ad7af79ce3cbc75015a28d4a8725786642 100644 (file)
@@ -974,6 +974,41 @@ namespace GridTools
   get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);
 
 
+  /**
+   * 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
index 13b2de6aa5f0b4ca1aa27dc4282fd7ba23635bda..b003cea19086389663f19bd7534929fa0e209f30 100644 (file)
@@ -2873,6 +2873,255 @@ next_cell:
 
 
 
+
+  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
    *
index dbb8abb26ebaa808b076365c26f4411fe97d6753..3aaeccbfe0384e48dc810c1a8868b019665204af 100644 (file)
@@ -224,7 +224,20 @@ 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<typename 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)

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.