]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduce new functions in TriaAccessor 10560/head
authorPeter Munch <peterrmuench@gmail.com>
Fri, 19 Jun 2020 18:10:02 +0000 (20:10 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Tue, 23 Jun 2020 17:57:58 +0000 (19:57 +0200)
doc/news/changes/minor/20200623Munch [new file with mode: 0644]
include/deal.II/grid/tria_accessor.h
include/deal.II/grid/tria_accessor.templates.h
source/dofs/dof_tools_constraints.cc
source/grid/grid_tools_dof_handlers.cc
source/grid/tria.cc
source/grid/tria_description.cc

diff --git a/doc/news/changes/minor/20200623Munch b/doc/news/changes/minor/20200623Munch
new file mode 100644 (file)
index 0000000..0e93ee2
--- /dev/null
@@ -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.
+<br>
+(Peter Munch, 2020/06/23)
index b1185cf998be0ebae9099f85f721f350eb7c2a92..703cc2ed6b80c02b13e52664e8970f5c4e0bcba6 100644 (file)
@@ -1605,6 +1605,49 @@ public:
   is_translation_of(
     const TriaIterator<TriaAccessor<structdim, dim, spacedim>> &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<unsigned int, unsigned int>
+  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<unsigned int, unsigned int>
+  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<unsigned int, unsigned int>
+  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<unsigned int, unsigned int>
+  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<unsigned int, unsigned int>
+  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<TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
-             GeometryInfo<dim>::faces_per_cell>
+  boost::container::small_vector<
+    TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
+    GeometryInfo<dim>::faces_per_cell>
   face_iterators() const;
 
   /**
index e60670ad3ff7955c7bbbd448965c65ee66c177d5..f65005cb9050ad81771a6312726b187925a9fbb0 100644 (file)
@@ -29,6 +29,8 @@
 #include <deal.II/grid/tria_iterator.templates.h>
 #include <deal.II/grid/tria_levels.h>
 
+#include <boost/container/small_vector.hpp>
+
 #include <cmath>
 
 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<structdim, dim, spacedim>::vertex_index(
   const unsigned int corner) const
 {
-  AssertIndexRange(corner, GeometryInfo<structdim>::vertices_per_cell);
+  AssertIndexRange(corner, this->n_vertices());
 
   return dealii::internal::TriaAccessorImplementation::Implementation::
     vertex_index(*this, corner);
@@ -1077,7 +1079,7 @@ template <int structdim, int dim, int spacedim>
 inline unsigned int
 TriaAccessor<structdim, dim, spacedim>::line_index(const unsigned int i) const
 {
-  AssertIndexRange(i, GeometryInfo<structdim>::lines_per_cell);
+  AssertIndexRange(i, this->n_lines());
 
   return dealii::internal::TriaAccessorImplementation::Implementation::
     line_index(*this, i);
@@ -1150,7 +1152,7 @@ TriaAccessor<structdim, dim, spacedim>::line_orientation(
   const unsigned int line) const
 {
   Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
-  AssertIndexRange(line, GeometryInfo<structdim>::lines_per_cell);
+  AssertIndexRange(line, this->n_lines());
 
   return dealii::internal::TriaAccessorImplementation::Implementation::
     line_orientation(*this, line);
@@ -1205,7 +1207,7 @@ TriaAccessor<structdim, dim, spacedim>::set_line_orientation(
   const bool         value) const
 {
   Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
-  AssertIndexRange(line, GeometryInfo<structdim>::lines_per_cell);
+  AssertIndexRange(line, this->n_lines());
 
   dealii::internal::TriaAccessorImplementation::Implementation::
     set_line_orientation(*this, line, value);
@@ -1945,8 +1947,7 @@ TriaAccessor<structdim, dim, spacedim>::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<bool, GeometryInfo<structdim>::vertices_per_cell>
-    is_initial_guess_vertex;
+  std::vector<bool> 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<structdim, dim, spacedim>::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<structdim>::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<structdim, dim, spacedim>::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<structdim>::vertex_indices())
+  for (const unsigned int v : this->vertex_indices())
     if (center.distance(this->vertex(v)) >
         radius + 100. * std::numeric_limits<double>::epsilon())
       {
@@ -2074,10 +2075,8 @@ TriaAccessor<structdim, dim, spacedim>::minimum_vertex_distance() const
       case 3:
         {
           double min = std::numeric_limits<double>::max();
-          for (const unsigned int i : GeometryInfo<structdim>::vertex_indices())
-            for (unsigned int j = i + 1;
-                 j < GeometryInfo<structdim>::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<structdim, dim, spacedim>::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<structdim>::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<structdim, dim, spacedim>::is_translation_of(
 }
 
 
+
+template <int structdim, int dim, int spacedim>
+unsigned int
+TriaAccessor<structdim, dim, spacedim>::n_vertices() const
+{
+  return GeometryInfo<structdim>::vertices_per_cell;
+}
+
+
+
+template <int structdim, int dim, int spacedim>
+unsigned int
+TriaAccessor<structdim, dim, spacedim>::n_lines() const
+{
+  return GeometryInfo<structdim>::lines_per_cell;
+}
+
+
+
+template <int structdim, int dim, int spacedim>
+unsigned int
+TriaAccessor<structdim, dim, spacedim>::n_faces() const
+{
+  AssertDimension(structdim, dim);
+
+  return GeometryInfo<structdim>::faces_per_cell;
+}
+
+
+
+template <int structdim, int dim, int spacedim>
+std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+TriaAccessor<structdim, dim, spacedim>::vertex_indices() const
+{
+  return {0U, n_vertices()};
+}
+
+
+
+template <int structdim, int dim, int spacedim>
+std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+TriaAccessor<structdim, dim, spacedim>::line_indices() const
+{
+  return {0U, n_lines()};
+}
+
+
+
+template <int structdim, int dim, int spacedim>
+std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+TriaAccessor<structdim, dim, spacedim>::face_indices() const
+{
+  return {0U, n_faces()};
+}
+
+
 /*----------------- Functions: TriaAccessor<0,dim,spacedim> -----------------*/
 
 template <int dim, int spacedim>
@@ -2974,6 +3029,42 @@ TriaAccessor<0, 1, spacedim>::used() const
   return tria->vertex_used(global_vertex_index);
 }
 
+
+
+template <int spacedim>
+unsigned int
+TriaAccessor<0, 1, spacedim>::n_vertices() const
+{
+  return 1;
+}
+
+
+
+template <int spacedim>
+unsigned int
+TriaAccessor<0, 1, spacedim>::n_lines() const
+{
+  return 0;
+}
+
+
+
+template <int spacedim>
+std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+TriaAccessor<0, 1, spacedim>::vertex_indices() const
+{
+  return {0U, n_vertices()};
+}
+
+
+
+template <int spacedim>
+std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+TriaAccessor<0, 1, spacedim>::line_indices() const
+{
+  return {0U, n_lines()};
+}
+
 /*------------------ Functions: CellAccessor<dim,spacedim> ------------------*/
 
 
@@ -3069,7 +3160,7 @@ inline unsigned int
 CellAccessor<dim, spacedim>::face_iterator_to_index(
   const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> &face) const
 {
-  for (const unsigned int face_n : GeometryInfo<dim>::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<dim, spacedim>::face_iterator_to_index(
 
 
 template <int dim, int spacedim>
-inline std::array<TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
-                  GeometryInfo<dim>::faces_per_cell>
+inline boost::container::small_vector<
+  TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
+  GeometryInfo<dim>::faces_per_cell>
 CellAccessor<dim, spacedim>::face_iterators() const
 {
-  std::array<TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
-             GeometryInfo<dim>::faces_per_cell>
-    face_iterators;
+  boost::container::small_vector<
+    TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
+    GeometryInfo<dim>::faces_per_cell>
+    face_iterators(this->n_faces());
 
-  for (unsigned int i : GeometryInfo<dim>::face_indices())
+  for (unsigned int i : this->face_indices())
     face_iterators[i] =
       dealii::internal::CellAccessorImplementation::get_face(*this, i);
 
@@ -3126,7 +3219,7 @@ template <int dim, int spacedim>
 inline int
 CellAccessor<dim, spacedim>::neighbor_index(const unsigned int i) const
 {
-  AssertIndexRange(i, GeometryInfo<dim>::faces_per_cell);
+  AssertIndexRange(i, this->n_faces());
   return this->tria->levels[this->present_level]
     ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
     .second;
@@ -3138,7 +3231,7 @@ template <int dim, int spacedim>
 inline int
 CellAccessor<dim, spacedim>::neighbor_level(const unsigned int i) const
 {
-  AssertIndexRange(i, GeometryInfo<dim>::faces_per_cell);
+  AssertIndexRange(i, this->n_faces());
   return this->tria->levels[this->present_level]
     ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
     .first;
@@ -3197,7 +3290,7 @@ CellAccessor<dim, spacedim>::flag_for_face_refinement(
   const RefinementCase<dim - 1> &face_refinement_case) const
 {
   Assert(dim > 1, ExcImpossibleInDim(dim));
-  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
+  AssertIndexRange(face_no, this->n_faces());
   AssertIndexRange(face_refinement_case,
                    RefinementCase<dim>::isotropic_refinement + 1);
 
@@ -3228,7 +3321,7 @@ CellAccessor<dim, spacedim>::flag_for_line_refinement(
   const unsigned int line_no) const
 {
   Assert(dim > 1, ExcImpossibleInDim(dim));
-  AssertIndexRange(line_no, GeometryInfo<dim>::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<std::uint8_t>(face(face_no)->refinement_case()))
     {
       case RefinementCase<3>::no_refinement:
index 708bcaaea73d3afa8fdb2229f7ee7c7a3d3c9bfd..72e63f0553ecdd7b65a9e05263e7bb5203cfd025 100644 (file)
@@ -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<dim>::face_indices())
+          for (const auto face_no : cell->face_indices())
             {
               const typename DoFHandlerType<dim, spacedim>::face_iterator face =
                 cell->face(face_no);
index 21c3c01faf7260d03bbaead7890b807ed0623625..2933623d74ab28e0f3540a5e73e02efad3bfdc03 100644 (file)
@@ -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<MeshType::dimension>::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<MeshType::dimension>::vertex_indices())
+        for (const auto v : cell->vertex_indices())
           if (locally_active_vertices_on_subdomain[cell->vertex_index(v)] ==
               true)
             {
index ce6d23b38f9e9de8ac7a93df0dff56d52cdb4ac1..f615882f23757da1cbe130d9d5cea6e8699450c7 100644 (file)
@@ -10972,16 +10972,12 @@ Triangulation<dim, spacedim>::create_triangulation(
             while (cell_info->id != cell->id().template to_binary<dim>())
               ++cell;
             if (dim == 3)
-              for (unsigned int quad = 0;
-                   quad < GeometryInfo<dim>::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<dim>::lines_per_cell;
-                   ++line)
+              for (const auto line : cell->line_indices())
                 cell->line(line)->set_manifold_id(
                   cell_info->manifold_line_ids[line]);
 
index bbcaab8b7b915d3d004fb9a8dd7b47c51dce7f2b..7dd871a6c515331cb8c7b8e688cd93f4e7e4a613 100644 (file)
@@ -161,7 +161,7 @@ namespace TriangulationDescription
           TriaIterator<CellAccessor<dim, spacedim>> &cell,
           std::vector<bool> &vertices_owned_by_locally_owned_cells) {
           // add local vertices
-          for (const auto v : GeometryInfo<dim>::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<CellAccessor<dim, spacedim>> &cell) {
-              for (const auto v : GeometryInfo<dim>::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<dim> cell_data(
-              GeometryInfo<dim>::vertices_per_cell);
+            dealii::CellData<dim> cell_data(cell->n_vertices());
             cell_data.material_id = cell->material_id();
             cell_data.manifold_id = cell->manifold_id();
-            for (const auto v : GeometryInfo<dim>::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<dim>::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<dim>::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<CellAccessor<dim, spacedim>> &cell) {
           if (cell->active())
-            for (const auto v : GeometryInfo<dim>::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<CellAccessor<dim, spacedim>> &cell) {
-              for (const auto v : GeometryInfo<dim>::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<dim>();
 
               // save boundary_ids of each face of this cell
-              for (const auto f : GeometryInfo<dim>::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<dim>::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<dim>::face_indices())
+                  for (const auto f : cell->face_indices())
                     cell_info.manifold_quad_ids[f] =
                       cell->quad(f)->manifold_id();
               }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.