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
+ 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
*
projects / dealii.git / commitdiff