]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add ReferenceCell::max_n_faces<dim>().
authorDavid Wells <drwells@email.unc.edu>
Mon, 30 Dec 2024 17:49:05 +0000 (12:49 -0500)
committerDavid Wells <drwells@email.unc.edu>
Tue, 31 Dec 2024 16:37:14 +0000 (11:37 -0500)
This is a more descriptive function than
GeometryInfo<dim>::faces_per_cell since in many cases we want to use the
maximum value (which happens to equal the hypercube value) regardless of
the reference cell type.

include/deal.II/grid/reference_cell.h
include/deal.II/grid/tria_accessor.h
include/deal.II/matrix_free/fe_evaluation_data.h
include/deal.II/matrix_free/mapping_info.templates.h
include/deal.II/matrix_free/matrix_free.h
include/deal.II/matrix_free/matrix_free.templates.h
source/grid/grid_in.cc
source/grid/tria.cc
source/grid/tria_accessor.cc

index df32f602bc220f79375407c3566d83687f4c0322..dca975c3adfb7204b301aadc16c572033d0d3dca 100644 (file)
@@ -326,6 +326,25 @@ public:
   unsigned int
   n_faces() const;
 
+  /**
+   * Return the maximum number of faces an object of dimension `structdim` can
+   * have. This is always the number of faces of a `structdim`-dimensional
+   * hypercube.
+   *
+   * @note The primary use case of this and the other maxima functions is to
+   * enable simple array indexing to per-face data by cell index and face
+   * number, e.g.,
+   *
+   * @code
+   *     cell->index() * ReferenceCell::max_n_faces<dim>() + face_no;
+   * @endcode
+   *
+   * is a unique index to the `face_no`th face of the current cell.
+   */
+  template <int structdim>
+  static constexpr unsigned int
+  max_n_faces();
+
   /**
    * Return an object that can be thought of as an array containing all
    * indices from zero to n_faces().
@@ -1741,6 +1760,15 @@ ReferenceCell::n_faces() const
 
 
 
+template <int structdim>
+inline constexpr unsigned int
+ReferenceCell::max_n_faces()
+{
+  return GeometryInfo<structdim>::faces_per_cell;
+}
+
+
+
 inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
 ReferenceCell::face_indices() const
 {
index 5eb6e412b0763f9e2f5c90ddd9c5611ff0fcd30b..4808a109c7f213d48e88954a9014e300b916b07a 100644 (file)
@@ -3240,7 +3240,7 @@ public:
    */
   boost::container::small_vector<
     TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
-    GeometryInfo<dim>::faces_per_cell>
+    ReferenceCell::max_n_faces<dim>()>
   face_iterators() const;
 
   /**
@@ -4912,9 +4912,9 @@ namespace internal
       inline static unsigned int
       quad_index(const TriaAccessor<3, 3, 3> &accessor, const unsigned int i)
       {
-        constexpr unsigned int max_faces_per_cell = 6;
         return accessor.tria->levels[accessor.present_level]
-          ->cells.cells[accessor.present_index * max_faces_per_cell + i];
+          ->cells
+          .cells[accessor.present_index * ReferenceCell::max_n_faces<3>() + i];
       }
 
 
@@ -4953,7 +4953,7 @@ namespace internal
         if (dim != 1)
           accessor.tria->levels[accessor.present_level]
             ->face_orientations.set_combined_orientation(
-              accessor.present_index * GeometryInfo<dim>::faces_per_cell + face,
+              accessor.present_index * ReferenceCell::max_n_faces<dim>() + face,
               combined_orientation);
       }
 
@@ -4967,9 +4967,9 @@ namespace internal
       vertex_index(const TriaAccessor<1, dim, spacedim> &accessor,
                    const unsigned int                    corner)
       {
-        constexpr unsigned int max_faces_per_cell = 2;
         return accessor.objects()
-          .cells[accessor.present_index * max_faces_per_cell + corner];
+          .cells[accessor.present_index * ReferenceCell::max_n_faces<1>() +
+                 corner];
       }
 
 
@@ -5059,7 +5059,7 @@ namespace internal
                   cell.get_triangulation()
                     .levels[cell.level()]
                     ->face_orientations.get_combined_orientation(
-                      cell.index() * GeometryInfo<3>::faces_per_cell + f);
+                      cell.index() * ReferenceCell::max_n_faces<dim>() + f);
 
                 // It might seem superfluous to spell out the four indices
                 // that get later consumed by a for loop over these four
@@ -5083,7 +5083,7 @@ namespace internal
                   cell.get_triangulation()
                     .levels[cell.level()]
                     ->face_orientations.get_combined_orientation(
-                      cell.index() * GeometryInfo<3>::faces_per_cell + f);
+                      cell.index() * ReferenceCell::max_n_faces<dim>() + f);
                 const std::array<unsigned int, 2> my_indices{
                   {ref_cell.standard_to_real_face_line(0, f, orientation),
                    ref_cell.standard_to_real_face_line(1, f, orientation)}};
@@ -5180,7 +5180,7 @@ namespace internal
                   cell.get_triangulation()
                     .levels[cell.level()]
                     ->face_orientations.get_combined_orientation(
-                      cell.index() * GeometryInfo<3>::faces_per_cell + f);
+                      cell.index() * ReferenceCell::max_n_faces<dim>() + f);
 
                 // It might seem superfluous to spell out the four indices and
                 // orientations that get later consumed by a for loop over
@@ -5215,7 +5215,7 @@ namespace internal
                   cell.get_triangulation()
                     .levels[cell.level()]
                     ->face_orientations.get_combined_orientation(
-                      cell.index() * GeometryInfo<3>::faces_per_cell + f);
+                      cell.index() * ReferenceCell::max_n_faces<3>() + f);
                 const std::array<unsigned int, 2> my_indices{
                   {ref_cell.standard_to_real_face_line(0, f, orientation),
                    ref_cell.standard_to_real_face_line(1, f, orientation)}};
@@ -5403,9 +5403,8 @@ TriaAccessor<structdim, dim, spacedim>::line_index(const unsigned int i) const
 
   if constexpr (structdim == 2)
     {
-      constexpr unsigned int max_faces_per_cell = 4;
       return this->objects()
-        .cells[this->present_index * max_faces_per_cell + i];
+        .cells[this->present_index * ReferenceCell::max_n_faces<2>() + i];
     }
   else if constexpr (structdim == 3)
     {
@@ -5472,13 +5471,12 @@ TriaAccessor<structdim, dim, spacedim>::combined_face_orientation(
       else
         return this->tria->levels[this->present_level]
           ->face_orientations.get_orientation(
-            this->present_index * GeometryInfo<structdim>::faces_per_cell +
-            face);
+            this->present_index * ReferenceCell::max_n_faces<dim>() + face);
     }
   else
     return this->tria->levels[this->present_level]
       ->face_orientations.get_combined_orientation(
-        this->present_index * GeometryInfo<structdim>::faces_per_cell + face);
+        this->present_index * ReferenceCell::max_n_faces<dim>() + face);
 }
 
 
@@ -5505,7 +5503,7 @@ TriaAccessor<structdim, dim, spacedim>::face_orientation(
   else
     return this->tria->levels[this->present_level]
       ->face_orientations.get_orientation(
-        this->present_index * GeometryInfo<structdim>::faces_per_cell + face);
+        this->present_index * ReferenceCell::max_n_faces<structdim>() + face);
 }
 
 
@@ -5525,7 +5523,7 @@ TriaAccessor<structdim, dim, spacedim>::face_flip(const unsigned int face) const
 
   if constexpr (structdim == 3)
     return this->tria->levels[this->present_level]->face_orientations.get_flip(
-      this->present_index * GeometryInfo<3>::faces_per_cell + face);
+      this->present_index * ReferenceCell::max_n_faces<structdim>() + face);
   else
     // In 1d and 2d, face_flip is always false as faces can only be
     // 'flipped' in 3d.
@@ -5550,7 +5548,7 @@ TriaAccessor<structdim, dim, spacedim>::face_rotation(
   if constexpr (structdim == 3)
     return this->tria->levels[this->present_level]
       ->face_orientations.get_rotation(
-        this->present_index * GeometryInfo<3>::faces_per_cell + face);
+        this->present_index * ReferenceCell::max_n_faces<structdim>() + face);
   else
     // In 1d and 2d, face_rotation is always false as faces can only be
     // 'rotated' in 3d.
@@ -7686,12 +7684,12 @@ CellAccessor<dim, spacedim>::face_iterator_to_index(
 template <int dim, int spacedim>
 inline boost::container::small_vector<
   TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
-  GeometryInfo<dim>::faces_per_cell>
+  ReferenceCell::max_n_faces<dim>()>
 CellAccessor<dim, spacedim>::face_iterators() const
 {
   boost::container::small_vector<
     TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
-    GeometryInfo<dim>::faces_per_cell>
+    ReferenceCell::max_n_faces<dim>()>
     face_iterators(this->n_faces());
 
   for (const unsigned int i : this->face_indices())
@@ -7732,7 +7730,7 @@ CellAccessor<dim, spacedim>::neighbor_index(const unsigned int face_no) const
 {
   AssertIndexRange(face_no, this->n_faces());
   return this->tria->levels[this->present_level]
-    ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell +
+    ->neighbors[this->present_index * ReferenceCell::max_n_faces<dim>() +
                 face_no]
     .second;
 }
@@ -7745,7 +7743,7 @@ CellAccessor<dim, spacedim>::neighbor_level(const unsigned int face_no) const
 {
   AssertIndexRange(face_no, this->n_faces());
   return this->tria->levels[this->present_level]
-    ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell +
+    ->neighbors[this->present_index * ReferenceCell::max_n_faces<dim>() +
                 face_no]
     .first;
 }
index de010f39ab073b147b03628ec488b13d447e761d..f4e3d2d8063a509eb38f25162a5a65874647f5aa 100644 (file)
@@ -30,6 +30,8 @@
 #include <deal.II/base/tensor.h>
 #include <deal.II/base/vectorization.h>
 
+#include <deal.II/grid/reference_cell.h>
+
 #include <deal.II/matrix_free/dof_info.h>
 #include <deal.II/matrix_free/mapping_info_storage.h>
 #include <deal.II/matrix_free/shape_info.h>
@@ -1542,7 +1544,7 @@ FEEvaluationData<dim, Number, is_face>::get_current_cell_index() const
 {
   if (is_face && dof_access_index ==
                    internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
-    return cell * GeometryInfo<dim>::faces_per_cell + face_numbers[0];
+    return cell * ReferenceCell::max_n_faces<dim>() + face_numbers[0];
   else
     return cell;
 }
index 00eacaf94ec78d968103dd7de89467c4147313a8..9b3c04d3b645974f24a893f8851bd4164eb73594 100644 (file)
@@ -3099,10 +3099,10 @@ namespace internal
           // counting
           AssertDimension(cell_type.size(), cells.size() / n_lanes);
           face_data_by_cells[my_q].data_index_offsets.resize(
-            cell_type.size() * GeometryInfo<dim>::faces_per_cell);
+            cell_type.size() * ReferenceCell::max_n_faces<dim>());
           if (update_flags & update_quadrature_points)
             face_data_by_cells[my_q].quadrature_point_offsets.resize(
-              cell_type.size() * GeometryInfo<dim>::faces_per_cell);
+              cell_type.size() * ReferenceCell::max_n_faces<dim>());
           std::size_t storage_length = 0;
           for (unsigned int i = 0; i < cell_type.size(); ++i)
             for (const unsigned int face : GeometryInfo<dim>::face_indices())
@@ -3110,54 +3110,54 @@ namespace internal
                 if (faces_by_cells_type[i][face] <= affine)
                   {
                     face_data_by_cells[my_q].data_index_offsets
-                      [i * GeometryInfo<dim>::faces_per_cell + face] =
+                      [i * ReferenceCell::max_n_faces<dim>() + face] =
                       storage_length;
                     ++storage_length;
                   }
                 else
                   {
                     face_data_by_cells[my_q].data_index_offsets
-                      [i * GeometryInfo<dim>::faces_per_cell + face] =
+                      [i * ReferenceCell::max_n_faces<dim>() + face] =
                       storage_length;
                     storage_length +=
                       face_data_by_cells[my_q].descriptor[0].n_q_points;
                   }
                 if (update_flags & update_quadrature_points)
                   face_data_by_cells[my_q].quadrature_point_offsets
-                    [i * GeometryInfo<dim>::faces_per_cell + face] =
-                    (i * GeometryInfo<dim>::faces_per_cell + face) *
+                    [i * ReferenceCell::max_n_faces<dim>() + face] =
+                    (i * ReferenceCell::max_n_faces<dim>() + face) *
                     face_data_by_cells[my_q].descriptor[0].n_q_points;
               }
           face_data_by_cells[my_q].JxW_values.resize_fast(
-            storage_length * GeometryInfo<dim>::faces_per_cell);
+            storage_length * ReferenceCell::max_n_faces<dim>());
           face_data_by_cells[my_q].jacobians[0].resize_fast(
-            storage_length * GeometryInfo<dim>::faces_per_cell);
+            storage_length * ReferenceCell::max_n_faces<dim>());
           face_data_by_cells[my_q].jacobians[1].resize_fast(
-            storage_length * GeometryInfo<dim>::faces_per_cell);
+            storage_length * ReferenceCell::max_n_faces<dim>());
           if (update_flags & update_normal_vectors)
             face_data_by_cells[my_q].normal_vectors.resize_fast(
-              storage_length * GeometryInfo<dim>::faces_per_cell);
+              storage_length * ReferenceCell::max_n_faces<dim>());
           if (update_flags & update_normal_vectors &&
               update_flags & update_jacobians)
             face_data_by_cells[my_q].normals_times_jacobians[0].resize_fast(
-              storage_length * GeometryInfo<dim>::faces_per_cell);
+              storage_length * ReferenceCell::max_n_faces<dim>());
           if (update_flags & update_normal_vectors &&
               update_flags & update_jacobians)
             face_data_by_cells[my_q].normals_times_jacobians[1].resize_fast(
-              storage_length * GeometryInfo<dim>::faces_per_cell);
+              storage_length * ReferenceCell::max_n_faces<dim>());
           if (update_flags & update_jacobian_grads)
             {
               face_data_by_cells[my_q].jacobian_gradients[0].resize_fast(
-                storage_length * GeometryInfo<dim>::faces_per_cell);
+                storage_length * ReferenceCell::max_n_faces<dim>());
               face_data_by_cells[my_q]
                 .jacobian_gradients_non_inverse[0]
                 .resize_fast(storage_length *
-                             GeometryInfo<dim>::faces_per_cell);
+                             ReferenceCell::max_n_faces<dim>());
             }
 
           if (update_flags & update_quadrature_points)
             face_data_by_cells[my_q].quadrature_points.resize_fast(
-              cell_type.size() * GeometryInfo<dim>::faces_per_cell *
+              cell_type.size() * ReferenceCell::max_n_faces<dim>() *
               face_data_by_cells[my_q].descriptor[0].n_q_points);
         }
 
@@ -3196,7 +3196,7 @@ namespace internal
                 *fe_face_values_neigh[my_q][fe_index];
               const unsigned int offset =
                 face_data_by_cells[my_q]
-                  .data_index_offsets[cell * GeometryInfo<dim>::faces_per_cell +
+                  .data_index_offsets[cell * ReferenceCell::max_n_faces<dim>() +
                                       face];
 
               const GeometryType my_cell_type = faces_by_cells_type[cell][face];
@@ -3334,7 +3334,7 @@ namespace internal
                       for (unsigned int d = 0; d < dim; ++d)
                         face_data_by_cells[my_q].quadrature_points
                           [face_data_by_cells[my_q].quadrature_point_offsets
-                             [cell * GeometryInfo<dim>::faces_per_cell + face] +
+                             [cell * ReferenceCell::max_n_faces<dim>() + face] +
                            q][d][v] = fe_val.quadrature_point(q)[d];
                 }
               if (update_flags & update_normal_vectors &&
@@ -3372,7 +3372,7 @@ namespace internal
       memory += cell_type.capacity() * sizeof(GeometryType);
       memory += face_type.capacity() * sizeof(GeometryType);
       memory += faces_by_cells_type.capacity() *
-                GeometryInfo<dim>::faces_per_cell * sizeof(GeometryType);
+                ReferenceCell::max_n_faces<dim>() * sizeof(GeometryType);
       memory += sizeof(*this);
       return memory;
     }
@@ -3398,7 +3398,7 @@ namespace internal
       out << "    Faces by cells types:            ";
       task_info.print_memory_statistics(out,
                                         faces_by_cells_type.capacity() *
-                                          GeometryInfo<dim>::faces_per_cell *
+                                          ReferenceCell::max_n_faces<dim>() *
                                           sizeof(GeometryType));
 
       for (unsigned int j = 0; j < cell_data.size(); ++j)
index a68bc1be24b73c1e8a8927a568170aa931047bd8..4eecd5539a14f71efe924ebdaa66e4288a87f2cb 100644 (file)
@@ -2571,7 +2571,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::get_faces_by_cells_boundary_id(
   const unsigned int face_number) const
 {
   AssertIndexRange(cell_batch_index, n_cell_batches());
-  AssertIndexRange(face_number, GeometryInfo<dim>::faces_per_cell);
+  AssertIndexRange(face_number, ReferenceCell::max_n_faces<dim>());
   Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(),
          ExcNotInitialized());
   std::array<types::boundary_id, VectorizedArrayType::size()> result;
index 11c72fd1b6bd6ff9c97b31c036f8ad3f310e00f9..23affe7c9028a6a6ed3dec4bd0b3b1a9f4c19c58 100644 (file)
@@ -2062,7 +2062,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
         numbers::invalid_unsigned_int);
       face_info.cell_and_face_boundary_id.reinit(
         TableIndices<3>(task_info.cell_partition_data.back(),
-                        GeometryInfo<dim>::faces_per_cell,
+                        ReferenceCell::max_n_faces<dim>(),
                         VectorizedArrayType::size()),
         true);
       face_info.cell_and_face_boundary_id.fill(numbers::invalid_boundary_id);
index 01523b61401a3c2ec644c9989f21b1840a0fb1ee..909902b64ea0823f50ddd1ef87c45eaa03c92294 100644 (file)
@@ -3737,7 +3737,6 @@ namespace
         // the side sets so that we can convert a set of side set indices into
         // a single deal.II boundary or manifold id (and save the
         // correspondence).
-        constexpr auto max_faces_per_cell = GeometryInfo<dim>::faces_per_cell;
         std::map<std::size_t, std::vector<int>> face_side_sets;
         for (const int side_set_id : side_set_ids)
           {
@@ -3772,7 +3771,7 @@ namespace
                     const long        element_n = elements[side_n] - 1;
                     const long        face_n    = faces[side_n] - 1;
                     const std::size_t face_id =
-                      element_n * max_faces_per_cell + face_n;
+                      element_n * ReferenceCell::max_n_faces<dim>() + face_n;
                     face_side_sets[face_id].push_back(side_set_id);
                   }
               }
@@ -3818,9 +3817,11 @@ namespace
                        ExcInternalError());
               }
             // Record the b_or_m_id of the current face.
-            const unsigned int   local_face_n = face_id % max_faces_per_cell;
-            const CellData<dim> &cell = cells[face_id / max_faces_per_cell];
-            const ReferenceCell  cell_type =
+            const unsigned int local_face_n =
+              face_id % ReferenceCell::max_n_faces<dim>();
+            const CellData<dim> &cell =
+              cells[face_id / ReferenceCell::max_n_faces<dim>()];
+            const ReferenceCell cell_type =
               ReferenceCell::n_vertices_to_type(dim, cell.vertices.size());
             const unsigned int deal_face_n =
               cell_type.exodusii_face_to_deal_face(local_face_n);
index 02305f7aee66ae25e4a9094053df1c1021f2d864..6eacddc97c03d85b95f18ea6dc8eeff71ae0b1bd 100644 (file)
@@ -1827,8 +1827,6 @@ namespace
           }
       }
   }
-
-
 } // end of anonymous namespace
 
 
@@ -3253,7 +3251,7 @@ namespace internal
             for (unsigned int line = 0; line < n_lines; ++line)
               for (unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
                    ++i, ++j)
-                lines_0.cells[line * GeometryInfo<1>::faces_per_cell + j] =
+                lines_0.cells[line * ReferenceCell::max_n_faces<1>() + j] =
                   crs.col[i]; // set vertex indices
           }
 
@@ -3354,18 +3352,18 @@ namespace internal
                 {
                   // set neighbor if not at boundary
                   if (nei.col[i] != static_cast<unsigned int>(-1))
-                    level.neighbors[cell * GeometryInfo<dim>::faces_per_cell +
+                    level.neighbors[cell * ReferenceCell::max_n_faces<dim>() +
                                     j] = {0, nei.col[i]};
 
                   // set face indices
-                  cells_0.cells[cell * GeometryInfo<dim>::faces_per_cell + j] =
+                  cells_0.cells[cell * ReferenceCell::max_n_faces<dim>() + j] =
                     crs.col[i];
 
                   // set face orientation if needed
                   if (orientation_needed)
                     {
                       level.face_orientations.set_combined_orientation(
-                        cell * GeometryInfo<dim>::faces_per_cell + j,
+                        cell * ReferenceCell::max_n_faces<dim>() + j,
                         connectivity.entity_orientations(dim - 1)
                           .get_combined_orientation(i));
                     }
@@ -15976,7 +15974,7 @@ void Triangulation<dim, spacedim>::reset_cell_vertex_indices_cache()
 
                 const unsigned char combined_orientation =
                   levels[l]->face_orientations.get_combined_orientation(
-                    cell->index() * GeometryInfo<3>::faces_per_cell + face);
+                    cell->index() * ReferenceCell::max_n_faces<dim>() + face);
                 std::array<unsigned int, 4> vertex_order{
                   {ref_cell.standard_to_real_face_vertex(0,
                                                          face,
index 7a2bf6f9c4834fc106c6632ec862709419f6763f..704ae97b0453d1362f49c6fbcefff79c47cc2e9b 100644 (file)
@@ -2360,19 +2360,19 @@ CellAccessor<dim, spacedim>::set_neighbor(
   if (pointer.state() == IteratorState::valid)
     {
       this->tria->levels[this->present_level]
-        ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
+        ->neighbors[this->present_index * ReferenceCell::max_n_faces<dim>() + i]
         .first = pointer->present_level;
       this->tria->levels[this->present_level]
-        ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
+        ->neighbors[this->present_index * ReferenceCell::max_n_faces<dim>() + i]
         .second = pointer->present_index;
     }
   else
     {
       this->tria->levels[this->present_level]
-        ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
+        ->neighbors[this->present_index * ReferenceCell::max_n_faces<dim>() + i]
         .first = -1;
       this->tria->levels[this->present_level]
-        ->neighbors[this->present_index * GeometryInfo<dim>::faces_per_cell + i]
+        ->neighbors[this->present_index * ReferenceCell::max_n_faces<dim>() + i]
         .second = -1;
     }
 }

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.