From: Spencer Patty Date: Wed, 3 Feb 2016 19:21:50 +0000 (-0600) Subject: add the get_dof_to_support_patch_map() function to grid_tools and instantiations. X-Git-Tag: v8.4.0-rc2~19^2~9 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8667f82b2ad9541736d5b2369153ea672346c5de;p=dealii.git add the get_dof_to_support_patch_map() function to grid_tools and instantiations. --- diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 268c96f378..1daf92ad7a 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -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 or + * parallel:distributed::Triangulation 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 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 or hp::DoFHandler. + * @param[in] dof_handler The Container which could be built on a Triangulation + * or a parallel:distributed::Triangulation 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 + std::map< types::global_dof_index,std::vector > + get_dof_to_support_patch_map(Container& dof_handler); + + /*@}*/ /** * @name Lower-dimensional meshes for parts of higher-dimensional meshes diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 13b2de6aa5..b003cea190 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2873,6 +2873,255 @@ next_cell: + + template + std::map< types::global_dof_index,std::vector > + 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 > dof_to_set_of_cells_map; + + std::vector local_dof_index; + std::vector local_face_dof_index; + std::vector 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 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::lines_per_cell; ++l) + if (cell->line(l)->has_children()) + for (unsigned int c=0; cline(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::faces_per_cell; ++f) + { + if (cell->face(f)->has_children()) + { + for (unsigned int c=0; cface(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 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; cneighbor(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::lines_per_cell; ++l) + { + if (cell->line(l)->has_children()) + { + for (unsigned int c=0; cline(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; iline(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; in_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; iline(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 > dof_to_cell_patches; + + typename std::map >::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 * diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index dbb8abb26e..3aaeccbfe0 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -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::active_cell_iterator> > + get_dof_to_support_patch_map > + (Container &dof_handler); + + + \} +#endif + } for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; Container : TRIANGULATION_AND_DOFHANDLER_TEMPLATES)