- template <int dim, template <int, int> class MeshType, int spacedim>
- DEAL_II_CXX20_REQUIRES(
- (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
-#ifndef _MSC_VER
- std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
-#else
- std::vector<
- typename dealii::internal::
- ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
-#endif
- find_cells_adjacent_to_vertex(const MeshType<dim, spacedim> &mesh,
- const unsigned int vertex)
- {
- // make sure that the given vertex is
- // an active vertex of the underlying
- // triangulation
- AssertIndexRange(vertex, mesh.get_triangulation().n_vertices());
- Assert(mesh.get_triangulation().get_used_vertices()[vertex],
- ExcVertexNotUsed(vertex));
-
- // use a set instead of a vector
- // to ensure that cells are inserted only
- // once
- std::set<typename dealii::internal::
- ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>
- adjacent_cells;
-
- // go through all active cells and look if the vertex is part of that cell
- //
- // in 1d, this is all we need to care about. in 2d/3d we also need to worry
- // that the vertex might be a hanging node on a face or edge of a cell; in
- // this case, we would want to add those cells as well on whose faces the
- // vertex is located but for which it is not a vertex itself.
- //
- // getting this right is a lot simpler in 2d than in 3d. in 2d, a hanging
- // node can only be in the middle of a face and we can query the neighboring
- // cell from the current cell. on the other hand, in 3d a hanging node
- // vertex can also be on an edge but there can be many other cells on
- // this edge and we can not access them from the cell we are currently
- // on.
- //
- // so, in the 3d case, if we run the algorithm as in 2d, we catch all
- // those cells for which the vertex we seek is on a *subface*, but we
- // miss the case of cells for which the vertex we seek is on a
- // sub-edge for which there is no corresponding sub-face (because the
- // immediate neighbor behind this face is not refined), see for example
- // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we
- // haven't yet found the vertex for the current cell we also need to
- // look at the mid-points of edges
- //
- // as a final note, deciding whether a neighbor is actually coarser is
- // simple in the case of isotropic refinement (we just need to look at
- // the level of the current and the neighboring cell). however, this
- // isn't so simple if we have used anisotropic refinement since then
- // the level of a cell is not indicative of whether it is coarser or
- // not than the current cell. ultimately, we want to add all cells on
- // which the vertex is, independent of whether they are coarser or
- // finer and so in the 2d case below we simply add *any* *active* neighbor.
- // in the worst case, we add cells multiple times to the adjacent_cells
- // list, but std::set throws out those cells already entered
- for (const auto &cell : mesh.active_cell_iterators())
- {
- for (const unsigned int v : cell->vertex_indices())
- if (cell->vertex_index(v) == vertex)
- {
- // OK, we found a cell that contains
- // the given vertex. We add it
- // to the list.
- adjacent_cells.insert(cell);
-
- // as explained above, in 2+d we need to check whether
- // this vertex is on a face behind which there is a
- // (possibly) coarser neighbor. if this is the case,
- // then we need to also add this neighbor
- if (dim >= 2)
- for (const auto face :
- cell->reference_cell().faces_for_given_vertex(v))
- if (!cell->at_boundary(face) &&
- cell->neighbor(face)->is_active())
- {
- // there is a (possibly) coarser cell behind a
- // face to which the vertex belongs. the
- // vertex we are looking at is then either a
- // vertex of that coarser neighbor, or it is a
- // hanging node on one of the faces of that
- // cell. in either case, it is adjacent to the
- // vertex, so add it to the list as well (if
- // the cell was already in the list then the
- // std::set makes sure that we get it only
- // once)
- adjacent_cells.insert(cell->neighbor(face));
- }
-
- // in any case, we have found a cell, so go to the next cell
- goto next_cell;
- }
-
- // in 3d also loop over the edges
- if (dim >= 3)
- {
- for (unsigned int e = 0; e < cell->n_lines(); ++e)
- if (cell->line(e)->has_children())
- // the only place where this vertex could have been
- // hiding is on the mid-edge point of the edge we
- // are looking at
- if (cell->line(e)->child(0)->vertex_index(1) == vertex)
- {
- adjacent_cells.insert(cell);
-
- // jump out of this tangle of nested loops
- goto next_cell;
- }
- }
-
- // in more than 3d we would probably have to do the same as
- // above also for even lower-dimensional objects
- Assert(dim <= 3, ExcNotImplemented());
-
- // move on to the next cell if we have found the
- // vertex on the current one
- next_cell:;
- }
-
- // if this was an active vertex then there needs to have been
- // at least one cell to which it is adjacent!
- Assert(adjacent_cells.size() > 0, ExcInternalError());
-
- // return the result as a vector, rather than the set we built above
- return std::vector<
- typename dealii::internal::
- ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type>(
- adjacent_cells.begin(), adjacent_cells.end());
- }
-
-
-
template <int dim, int spacedim>
std::vector<std::vector<Tensor<1, spacedim>>>
vertex_to_cell_centers_directions(
projects / dealii.git / commitdiff