From: Peter Munch Date: Fri, 19 Jun 2020 18:10:02 +0000 (+0200) Subject: Introduce new functions in TriaAccessor X-Git-Tag: v9.3.0-rc1~1391^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10560%2Fhead;p=dealii.git Introduce new functions in TriaAccessor --- diff --git a/doc/news/changes/minor/20200623Munch b/doc/news/changes/minor/20200623Munch new file mode 100644 index 0000000000..0e93ee2f0b --- /dev/null +++ b/doc/news/changes/minor/20200623Munch @@ -0,0 +1,7 @@ +New: The class TriaAccessor provides now the capability to query +the number of vertices, lines, and faces (with n_vertices(), +n_lines(), n_faces(), vertex_indices(), +line_indices(), face_indices()). The new methods can be used as an +alternative to the methods in GeometryInfo. +
+(Peter Munch, 2020/06/23) diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index b1185cf998..703cc2ed6b 100644 --- a/include/deal.II/grid/tria_accessor.h +++ b/include/deal.II/grid/tria_accessor.h @@ -1605,6 +1605,49 @@ public: is_translation_of( const TriaIterator> &o) const; + /** + * Number of vertices. + */ + inline unsigned int + n_vertices() const; + + /** + * Number of lines. + */ + inline unsigned int + n_lines() const; + + /** + * Number of faces. + * + * @note Only implemented for cells (dim==spacedim). + */ + inline unsigned int + n_faces() const; + + /** + * Return an object that can be thought of as an array containing all indices + * from zero to n_vertices(). + */ + inline std_cxx20::ranges::iota_view + vertex_indices() const; + + /** + * Return an object that can be thought of as an array containing all indices + * from zero to n_lines(). + */ + inline std_cxx20::ranges::iota_view + line_indices() const; + + /** + * Return an object that can be thought of as an array containing all indices + * from zero to n_faces(). + * + * @note Only implemented for cells (dim==spacedim). + */ + inline std_cxx20::ranges::iota_view + face_indices() const; + /** * @} */ @@ -2628,6 +2671,32 @@ public: bool used() const; + /** + * Number of vertices. + */ + inline unsigned int + n_vertices() const; + + /** + * Number of lines. + */ + inline unsigned int + n_lines() const; + + /** + * Return an object that can be thought of as an array containing all indices + * from zero to n_vertices(). + */ + inline std_cxx20::ranges::iota_view + vertex_indices() const; + + /** + * Return an object that can be thought of as an array containing all indices + * from zero to n_lines(). + */ + inline std_cxx20::ranges::iota_view + line_indices() const; + protected: /** * Pointer to the triangulation we operate on. @@ -2766,8 +2835,9 @@ public: /** * Return an array of iterators to all faces of this cell. */ - std::array>, - GeometryInfo::faces_per_cell> + boost::container::small_vector< + TriaIterator>, + GeometryInfo::faces_per_cell> face_iterators() const; /** diff --git a/include/deal.II/grid/tria_accessor.templates.h b/include/deal.II/grid/tria_accessor.templates.h index e60670ad3f..f65005cb90 100644 --- a/include/deal.II/grid/tria_accessor.templates.h +++ b/include/deal.II/grid/tria_accessor.templates.h @@ -29,6 +29,8 @@ #include #include +#include + #include DEAL_II_NAMESPACE_OPEN @@ -660,7 +662,7 @@ namespace internal inline static bool face_flip(const TriaAccessor<3, 3, 3> &accessor, const unsigned int face) { - AssertIndexRange(face, GeometryInfo<3>::faces_per_cell); + AssertIndexRange(face, accessor.n_faces()); Assert(accessor.present_index * GeometryInfo<3>::faces_per_cell + face < accessor.tria->levels[accessor.present_level] ->face_orientations.size(), @@ -696,7 +698,7 @@ namespace internal face_rotation(const TriaAccessor<3, 3, 3> &accessor, const unsigned int face) { - AssertIndexRange(face, GeometryInfo<3>::faces_per_cell); + AssertIndexRange(face, accessor.n_faces()); Assert(accessor.present_index * GeometryInfo<3>::faces_per_cell + face < accessor.tria->levels[accessor.present_level] ->face_orientations.size(), @@ -750,7 +752,7 @@ namespace internal const unsigned int line) { Assert(accessor.used(), TriaAccessorExceptions::ExcCellNotUsed()); - AssertIndexRange(line, GeometryInfo<3>::lines_per_cell); + AssertIndexRange(line, accessor.n_lines()); const auto pair = GeometryInfo<3>::standard_hex_line_to_quad_line_index(line); @@ -821,7 +823,7 @@ namespace internal const bool value) { Assert(accessor.used(), TriaAccessorExceptions::ExcCellNotUsed()); - AssertIndexRange(face, GeometryInfo<3>::faces_per_cell); + AssertIndexRange(face, accessor.n_faces()); Assert(accessor.present_index * GeometryInfo<3>::faces_per_cell + face < accessor.tria->levels[accessor.present_level] ->face_orientations.size(), @@ -853,7 +855,7 @@ namespace internal const unsigned int face, const bool value) { - AssertIndexRange(face, GeometryInfo<3>::faces_per_cell); + AssertIndexRange(face, accessor.n_faces()); Assert(accessor.present_index * GeometryInfo<3>::faces_per_cell + face < accessor.tria->levels[accessor.present_level] ->face_orientations.size(), @@ -886,7 +888,7 @@ namespace internal const unsigned int face, const bool value) { - AssertIndexRange(face, GeometryInfo<3>::faces_per_cell); + AssertIndexRange(face, accessor.n_faces()); Assert(accessor.present_index * GeometryInfo<3>::faces_per_cell + face < accessor.tria->levels[accessor.present_level] ->face_orientations.size(), @@ -931,14 +933,14 @@ namespace internal const bool value) { Assert(accessor.used(), TriaAccessorExceptions::ExcCellNotUsed()); - AssertIndexRange(line, GeometryInfo<3>::lines_per_face); - Assert(accessor.present_index * GeometryInfo<3>::lines_per_face + line < + AssertIndexRange(line, accessor.n_lines()); + Assert(accessor.present_index * GeometryInfo<2>::lines_per_cell + line < accessor.tria->faces->quads_line_orientations.size(), ExcInternalError()); // quads as part of 3d hexes can have non-standard orientation accessor.tria->faces->quads_line_orientations - [accessor.present_index * GeometryInfo<3>::lines_per_face + line] = + [accessor.present_index * GeometryInfo<2>::lines_per_cell + line] = value; } @@ -1044,7 +1046,7 @@ inline unsigned int TriaAccessor::vertex_index( const unsigned int corner) const { - AssertIndexRange(corner, GeometryInfo::vertices_per_cell); + AssertIndexRange(corner, this->n_vertices()); return dealii::internal::TriaAccessorImplementation::Implementation:: vertex_index(*this, corner); @@ -1077,7 +1079,7 @@ template inline unsigned int TriaAccessor::line_index(const unsigned int i) const { - AssertIndexRange(i, GeometryInfo::lines_per_cell); + AssertIndexRange(i, this->n_lines()); return dealii::internal::TriaAccessorImplementation::Implementation:: line_index(*this, i); @@ -1150,7 +1152,7 @@ TriaAccessor::line_orientation( const unsigned int line) const { Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); - AssertIndexRange(line, GeometryInfo::lines_per_cell); + AssertIndexRange(line, this->n_lines()); return dealii::internal::TriaAccessorImplementation::Implementation:: line_orientation(*this, line); @@ -1205,7 +1207,7 @@ TriaAccessor::set_line_orientation( const bool value) const { Assert(used(), TriaAccessorExceptions::ExcCellNotUsed()); - AssertIndexRange(line, GeometryInfo::lines_per_cell); + AssertIndexRange(line, this->n_lines()); dealii::internal::TriaAccessorImplementation::Implementation:: set_line_orientation(*this, line, value); @@ -1945,8 +1947,7 @@ TriaAccessor::enclosing_ball() const // two vertices corresponding to the largest diagonal which is being used // to construct the initial ball. // We employ this mask to skip these two vertices while enlarging the ball. - std::array::vertices_per_cell> - is_initial_guess_vertex; + std::vector is_initial_guess_vertex(this->n_vertices()); // First let all the vertices be outside std::fill(is_initial_guess_vertex.begin(), @@ -2026,7 +2027,7 @@ TriaAccessor::enclosing_ball() const // For each vertex that is found to be geometrically outside the ball // enlarge the ball so that the new ball contains both the previous ball // and the given vertex. - for (const unsigned int v : GeometryInfo::vertex_indices()) + for (const unsigned int v : this->vertex_indices()) if (!is_initial_guess_vertex[v]) { const double distance = center.distance(this->vertex(v)); @@ -2048,7 +2049,7 @@ TriaAccessor::enclosing_ball() const // Set all_vertices_within_ball false if any of the vertices of the object // are geometrically outside the ball - for (const unsigned int v : GeometryInfo::vertex_indices()) + for (const unsigned int v : this->vertex_indices()) if (center.distance(this->vertex(v)) > radius + 100. * std::numeric_limits::epsilon()) { @@ -2074,10 +2075,8 @@ TriaAccessor::minimum_vertex_distance() const case 3: { double min = std::numeric_limits::max(); - for (const unsigned int i : GeometryInfo::vertex_indices()) - for (unsigned int j = i + 1; - j < GeometryInfo::vertices_per_cell; - ++j) + for (const unsigned int i : this->vertex_indices()) + for (unsigned int j = i + 1; j < this->n_vertices(); ++j) min = std::min(min, (this->vertex(i) - this->vertex(j)) * (this->vertex(i) - this->vertex(j))); @@ -2113,7 +2112,7 @@ TriaAccessor::is_translation_of( bool is_translation = true; const Tensor<1, spacedim> dist = o->vertex(0) - this->vertex(0); const double tol_square = 1e-24 * dist.norm_square(); - for (unsigned int i = 1; i < GeometryInfo::vertices_per_cell; ++i) + for (unsigned int i = 1; i < this->n_vertices(); ++i) { const Tensor<1, spacedim> dist_new = (o->vertex(i) - this->vertex(i)) - dist; @@ -2127,6 +2126,62 @@ TriaAccessor::is_translation_of( } + +template +unsigned int +TriaAccessor::n_vertices() const +{ + return GeometryInfo::vertices_per_cell; +} + + + +template +unsigned int +TriaAccessor::n_lines() const +{ + return GeometryInfo::lines_per_cell; +} + + + +template +unsigned int +TriaAccessor::n_faces() const +{ + AssertDimension(structdim, dim); + + return GeometryInfo::faces_per_cell; +} + + + +template +std_cxx20::ranges::iota_view +TriaAccessor::vertex_indices() const +{ + return {0U, n_vertices()}; +} + + + +template +std_cxx20::ranges::iota_view +TriaAccessor::line_indices() const +{ + return {0U, n_lines()}; +} + + + +template +std_cxx20::ranges::iota_view +TriaAccessor::face_indices() const +{ + return {0U, n_faces()}; +} + + /*----------------- Functions: TriaAccessor<0,dim,spacedim> -----------------*/ template @@ -2974,6 +3029,42 @@ TriaAccessor<0, 1, spacedim>::used() const return tria->vertex_used(global_vertex_index); } + + +template +unsigned int +TriaAccessor<0, 1, spacedim>::n_vertices() const +{ + return 1; +} + + + +template +unsigned int +TriaAccessor<0, 1, spacedim>::n_lines() const +{ + return 0; +} + + + +template +std_cxx20::ranges::iota_view +TriaAccessor<0, 1, spacedim>::vertex_indices() const +{ + return {0U, n_vertices()}; +} + + + +template +std_cxx20::ranges::iota_view +TriaAccessor<0, 1, spacedim>::line_indices() const +{ + return {0U, n_lines()}; +} + /*------------------ Functions: CellAccessor ------------------*/ @@ -3069,7 +3160,7 @@ inline unsigned int CellAccessor::face_iterator_to_index( const TriaIterator> &face) const { - for (const unsigned int face_n : GeometryInfo::face_indices()) + for (const unsigned int face_n : this->face_indices()) if (this->face(face_n) == face) return face_n; @@ -3081,15 +3172,17 @@ CellAccessor::face_iterator_to_index( template -inline std::array>, - GeometryInfo::faces_per_cell> +inline boost::container::small_vector< + TriaIterator>, + GeometryInfo::faces_per_cell> CellAccessor::face_iterators() const { - std::array>, - GeometryInfo::faces_per_cell> - face_iterators; + boost::container::small_vector< + TriaIterator>, + GeometryInfo::faces_per_cell> + face_iterators(this->n_faces()); - for (unsigned int i : GeometryInfo::face_indices()) + for (unsigned int i : this->face_indices()) face_iterators[i] = dealii::internal::CellAccessorImplementation::get_face(*this, i); @@ -3126,7 +3219,7 @@ template inline int CellAccessor::neighbor_index(const unsigned int i) const { - AssertIndexRange(i, GeometryInfo::faces_per_cell); + AssertIndexRange(i, this->n_faces()); return this->tria->levels[this->present_level] ->neighbors[this->present_index * GeometryInfo::faces_per_cell + i] .second; @@ -3138,7 +3231,7 @@ template inline int CellAccessor::neighbor_level(const unsigned int i) const { - AssertIndexRange(i, GeometryInfo::faces_per_cell); + AssertIndexRange(i, this->n_faces()); return this->tria->levels[this->present_level] ->neighbors[this->present_index * GeometryInfo::faces_per_cell + i] .first; @@ -3197,7 +3290,7 @@ CellAccessor::flag_for_face_refinement( const RefinementCase &face_refinement_case) const { Assert(dim > 1, ExcImpossibleInDim(dim)); - AssertIndexRange(face_no, GeometryInfo::faces_per_cell); + AssertIndexRange(face_no, this->n_faces()); AssertIndexRange(face_refinement_case, RefinementCase::isotropic_refinement + 1); @@ -3228,7 +3321,7 @@ CellAccessor::flag_for_line_refinement( const unsigned int line_no) const { Assert(dim > 1, ExcImpossibleInDim(dim)); - AssertIndexRange(line_no, GeometryInfo::lines_per_cell); + AssertIndexRange(line_no, this->n_lines()); // the new refinement case is a combination // of the minimum required one for the given @@ -3275,7 +3368,7 @@ inline dealii::internal::SubfaceCase<2> CellAccessor<2>::subface_case(const unsigned int face_no) const { Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive()); - AssertIndexRange(face_no, GeometryInfo<2>::faces_per_cell); + AssertIndexRange(face_no, this->n_faces()); return ((face(face_no)->has_children()) ? dealii::internal::SubfaceCase<2>::case_x : dealii::internal::SubfaceCase<2>::case_none); @@ -3286,7 +3379,7 @@ inline dealii::internal::SubfaceCase<2> CellAccessor<2, 3>::subface_case(const unsigned int face_no) const { Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive()); - AssertIndexRange(face_no, GeometryInfo<2>::faces_per_cell); + AssertIndexRange(face_no, this->n_faces()); return ((face(face_no)->has_children()) ? dealii::internal::SubfaceCase<2>::case_x : dealii::internal::SubfaceCase<2>::case_none); @@ -3298,7 +3391,7 @@ inline dealii::internal::SubfaceCase<3> CellAccessor<3>::subface_case(const unsigned int face_no) const { Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive()); - AssertIndexRange(face_no, GeometryInfo<3>::faces_per_cell); + AssertIndexRange(face_no, this->n_faces()); switch (static_cast(face(face_no)->refinement_case())) { case RefinementCase<3>::no_refinement: diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index 708bcaaea7..72e63f0553 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -3384,7 +3384,7 @@ namespace DoFTools cell_dofs.resize(fe.dofs_per_cell); cell->get_dof_indices(cell_dofs); - for (const unsigned int face_no : GeometryInfo::face_indices()) + for (const auto face_no : cell->face_indices()) { const typename DoFHandlerType::face_iterator face = cell->face(face_no); diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index 21c3c01faf..2933623d74 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -757,8 +757,7 @@ namespace GridTools // the halo layer for (const auto &cell : mesh.active_cell_iterators()) if (predicate(cell)) // True predicate --> Part of subdomain - for (const unsigned int v : - GeometryInfo::vertex_indices()) + for (const auto v : cell->vertex_indices()) locally_active_vertices_on_subdomain[cell->vertex_index(v)] = true; // Find the cells that do not conform to the predicate @@ -766,8 +765,7 @@ namespace GridTools // These comprise the halo layer for (const auto &cell : mesh.active_cell_iterators()) if (!predicate(cell)) // False predicate --> Potential halo cell - for (const unsigned int v : - GeometryInfo::vertex_indices()) + for (const auto v : cell->vertex_indices()) if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] == true) { diff --git a/source/grid/tria.cc b/source/grid/tria.cc index ce6d23b38f..f615882f23 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -10972,16 +10972,12 @@ Triangulation::create_triangulation( while (cell_info->id != cell->id().template to_binary()) ++cell; if (dim == 3) - for (unsigned int quad = 0; - quad < GeometryInfo::quads_per_cell; - ++quad) + for (const auto quad : cell->face_indices()) cell->quad(quad)->set_manifold_id( cell_info->manifold_quad_ids[quad]); if (dim >= 2) - for (unsigned int line = 0; - line < GeometryInfo::lines_per_cell; - ++line) + for (const auto line : cell->line_indices()) cell->line(line)->set_manifold_id( cell_info->manifold_line_ids[line]); diff --git a/source/grid/tria_description.cc b/source/grid/tria_description.cc index bbcaab8b7b..7dd871a6c5 100644 --- a/source/grid/tria_description.cc +++ b/source/grid/tria_description.cc @@ -161,7 +161,7 @@ namespace TriangulationDescription TriaIterator> &cell, std::vector &vertices_owned_by_locally_owned_cells) { // add local vertices - for (const auto v : GeometryInfo::vertex_indices()) + for (const auto v : cell->vertex_indices()) { vertices_owned_by_locally_owned_cells[cell->vertex_index(v)] = true; @@ -223,9 +223,9 @@ namespace TriangulationDescription // helper function to determine if cell is locally relevant // (i.e. a cell which is connected to a vertex via a locally // owned cell) - auto is_locally_relevant_on_level = + const auto is_locally_relevant_on_level = [&](TriaIterator> &cell) { - for (const auto v : GeometryInfo::vertex_indices()) + for (const auto v : cell->vertex_indices()) if (vertices_owned_by_locally_owned_cells_on_level [cell->vertex_index(v)]) return true; @@ -249,16 +249,15 @@ namespace TriangulationDescription continue; // extract cell definition (with old numbering of vertices) - dealii::CellData cell_data( - GeometryInfo::vertices_per_cell); + dealii::CellData cell_data(cell->n_vertices()); cell_data.material_id = cell->material_id(); cell_data.manifold_id = cell->manifold_id(); - for (const auto v : GeometryInfo::vertex_indices()) + for (const auto v : cell->vertex_indices()) cell_data.vertices[v] = cell->vertex_index(v); construction_data.coarse_cells.push_back(cell_data); // save indices of each vertex of this cell - for (const auto v : GeometryInfo::vertex_indices()) + for (const auto v : cell->vertex_indices()) vertices_locally_relevant[cell->vertex_index(v)] = numbers::invalid_unsigned_int; @@ -278,7 +277,7 @@ namespace TriangulationDescription // c) correct vertices of cells (make them local) for (auto &cell : construction_data.coarse_cells) - for (const auto v : GeometryInfo::vertex_indices()) + for (unsigned int v = 0; v < cell.vertices.size(); ++v) cell.vertices[v] = vertices_locally_relevant[cell.vertices[v]]; } @@ -297,10 +296,10 @@ namespace TriangulationDescription // helper function to determine if cell is locally relevant // on active level - auto is_locally_relevant_on_active_level = + const auto is_locally_relevant_on_active_level = [&](TriaIterator> &cell) { if (cell->active()) - for (const auto v : GeometryInfo::vertex_indices()) + for (const auto v : cell->vertex_indices()) if (vertices_owned_by_locally_owned_active_cells [cell->vertex_index(v)]) return true; @@ -323,9 +322,9 @@ namespace TriangulationDescription // helper function to determine if cell is locally relevant // on level - auto is_locally_relevant_on_level = + const auto is_locally_relevant_on_level = [&](TriaIterator> &cell) { - for (const auto v : GeometryInfo::vertex_indices()) + for (const auto v : cell->vertex_indices()) if (vertices_owned_by_locally_owned_cells_on_level [cell->vertex_index(v)]) return true; @@ -345,7 +344,7 @@ namespace TriangulationDescription cell_info.id = cell->id().template to_binary(); // save boundary_ids of each face of this cell - for (const auto f : GeometryInfo::face_indices()) + for (const auto f : cell->face_indices()) { types::boundary_id boundary_ind = cell->face(f)->boundary_id(); @@ -360,15 +359,13 @@ namespace TriangulationDescription // ... of lines if (dim >= 2) - for (unsigned int line = 0; - line < GeometryInfo::lines_per_cell; - ++line) + for (const auto line : cell->line_indices()) cell_info.manifold_line_ids[line] = cell->line(line)->manifold_id(); // ... of quads if (dim == 3) - for (const auto f : GeometryInfo::face_indices()) + for (const auto f : cell->face_indices()) cell_info.manifold_quad_ids[f] = cell->quad(f)->manifold_id(); }