From: Spencer Patty 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 Container 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 + std::vector + get_cells_at_coarsest_common_level(const std::vector &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 element. + * + * @tparam Container In C++, the compiler can not determine the type + * of Container 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 + void + build_triangulation_from_patch ( + const std::vector &patch, + Triangulation &local_triangulation, + std::map::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 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 --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 + std::vector + get_cells_at_coarsest_common_level ( + const std::vector &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; ilevel() ); + std::set uniform_cells; + typename std::vector::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 (uniform_cells.begin(), + uniform_cells.end()); + } + + + + template + void build_triangulation_from_patch(const std::vector &patch, + Triangulation &local_triangulation, + std::map::active_cell_iterator, + typename Container::active_cell_iterator> &patch_to_global_tria_map) + + { + const std::vector uniform_cells = + get_cells_at_coarsest_common_level (patch); + // First it creates triangulation from the vector of "uniform_cells" + local_triangulation.clear(); + std::vector > vertices; + const unsigned int n_uniform_cells=uniform_cells.size(); + std::vector > cells(n_uniform_cells); + unsigned int k=0;// for enumerating cells + unsigned int i=0;// for enumerating vertices + typename std::vector::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::vertices_per_cell; ++j) + { + Point position=(*uniform_cell)->vertex (j); + repeat_vertex=false; + + for (unsigned int m=0; m::cell_iterator, + typename Container::cell_iterator> patch_to_global_tria_map_tmp; + for (typename Triangulation::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::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; iset_user_flag(); + break; + } + } + } + + if (refinement_necessary) + { + local_triangulation.execute_coarsening_and_refinement (); + + for (typename Triangulation::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::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 + 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 *