From 411c749d9b57e01582785da49ebda3376b969283 Mon Sep 17 00:00:00 2001 From: David Wells Date: Mon, 30 Dec 2024 12:49:05 -0500 Subject: [PATCH] Add ReferenceCell::max_n_faces(). This is a more descriptive function than GeometryInfo::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 | 28 +++++++++++++ include/deal.II/grid/tria_accessor.h | 42 +++++++++---------- .../deal.II/matrix_free/fe_evaluation_data.h | 4 +- .../matrix_free/mapping_info.templates.h | 38 ++++++++--------- include/deal.II/matrix_free/matrix_free.h | 2 +- .../matrix_free/matrix_free.templates.h | 2 +- source/grid/grid_in.cc | 11 ++--- source/grid/tria.cc | 12 +++--- source/grid/tria_accessor.cc | 8 ++-- 9 files changed, 87 insertions(+), 60 deletions(-) diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index df32f602bc..dca975c3ad 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -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() + face_no; + * @endcode + * + * is a unique index to the `face_no`th face of the current cell. + */ + template + 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 +inline constexpr unsigned int +ReferenceCell::max_n_faces() +{ + return GeometryInfo::faces_per_cell; +} + + + inline std_cxx20::ranges::iota_view ReferenceCell::face_indices() const { diff --git a/include/deal.II/grid/tria_accessor.h b/include/deal.II/grid/tria_accessor.h index 5eb6e412b0..4808a109c7 100644 --- a/include/deal.II/grid/tria_accessor.h +++ b/include/deal.II/grid/tria_accessor.h @@ -3240,7 +3240,7 @@ public: */ boost::container::small_vector< TriaIterator>, - GeometryInfo::faces_per_cell> + ReferenceCell::max_n_faces()> 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::faces_per_cell + face, + accessor.present_index * ReferenceCell::max_n_faces() + 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() + 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() + f); const std::array 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() + 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 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::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::combined_face_orientation( else return this->tria->levels[this->present_level] ->face_orientations.get_orientation( - this->present_index * GeometryInfo::faces_per_cell + - face); + this->present_index * ReferenceCell::max_n_faces() + face); } else return this->tria->levels[this->present_level] ->face_orientations.get_combined_orientation( - this->present_index * GeometryInfo::faces_per_cell + face); + this->present_index * ReferenceCell::max_n_faces() + face); } @@ -5505,7 +5503,7 @@ TriaAccessor::face_orientation( else return this->tria->levels[this->present_level] ->face_orientations.get_orientation( - this->present_index * GeometryInfo::faces_per_cell + face); + this->present_index * ReferenceCell::max_n_faces() + face); } @@ -5525,7 +5523,7 @@ TriaAccessor::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() + face); else // In 1d and 2d, face_flip is always false as faces can only be // 'flipped' in 3d. @@ -5550,7 +5548,7 @@ TriaAccessor::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() + face); else // In 1d and 2d, face_rotation is always false as faces can only be // 'rotated' in 3d. @@ -7686,12 +7684,12 @@ CellAccessor::face_iterator_to_index( template inline boost::container::small_vector< TriaIterator>, - GeometryInfo::faces_per_cell> + ReferenceCell::max_n_faces()> CellAccessor::face_iterators() const { boost::container::small_vector< TriaIterator>, - GeometryInfo::faces_per_cell> + ReferenceCell::max_n_faces()> face_iterators(this->n_faces()); for (const unsigned int i : this->face_indices()) @@ -7732,7 +7730,7 @@ CellAccessor::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::faces_per_cell + + ->neighbors[this->present_index * ReferenceCell::max_n_faces() + face_no] .second; } @@ -7745,7 +7743,7 @@ CellAccessor::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::faces_per_cell + + ->neighbors[this->present_index * ReferenceCell::max_n_faces() + face_no] .first; } diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h index de010f39ab..f4e3d2d806 100644 --- a/include/deal.II/matrix_free/fe_evaluation_data.h +++ b/include/deal.II/matrix_free/fe_evaluation_data.h @@ -30,6 +30,8 @@ #include #include +#include + #include #include #include @@ -1542,7 +1544,7 @@ FEEvaluationData::get_current_cell_index() const { if (is_face && dof_access_index == internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) - return cell * GeometryInfo::faces_per_cell + face_numbers[0]; + return cell * ReferenceCell::max_n_faces() + face_numbers[0]; else return cell; } diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 00eacaf94e..9b3c04d3b6 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -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::faces_per_cell); + cell_type.size() * ReferenceCell::max_n_faces()); if (update_flags & update_quadrature_points) face_data_by_cells[my_q].quadrature_point_offsets.resize( - cell_type.size() * GeometryInfo::faces_per_cell); + cell_type.size() * ReferenceCell::max_n_faces()); std::size_t storage_length = 0; for (unsigned int i = 0; i < cell_type.size(); ++i) for (const unsigned int face : GeometryInfo::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::faces_per_cell + face] = + [i * ReferenceCell::max_n_faces() + face] = storage_length; ++storage_length; } else { face_data_by_cells[my_q].data_index_offsets - [i * GeometryInfo::faces_per_cell + face] = + [i * ReferenceCell::max_n_faces() + 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::faces_per_cell + face] = - (i * GeometryInfo::faces_per_cell + face) * + [i * ReferenceCell::max_n_faces() + face] = + (i * ReferenceCell::max_n_faces() + 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::faces_per_cell); + storage_length * ReferenceCell::max_n_faces()); face_data_by_cells[my_q].jacobians[0].resize_fast( - storage_length * GeometryInfo::faces_per_cell); + storage_length * ReferenceCell::max_n_faces()); face_data_by_cells[my_q].jacobians[1].resize_fast( - storage_length * GeometryInfo::faces_per_cell); + storage_length * ReferenceCell::max_n_faces()); if (update_flags & update_normal_vectors) face_data_by_cells[my_q].normal_vectors.resize_fast( - storage_length * GeometryInfo::faces_per_cell); + storage_length * ReferenceCell::max_n_faces()); 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::faces_per_cell); + storage_length * ReferenceCell::max_n_faces()); 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::faces_per_cell); + storage_length * ReferenceCell::max_n_faces()); if (update_flags & update_jacobian_grads) { face_data_by_cells[my_q].jacobian_gradients[0].resize_fast( - storage_length * GeometryInfo::faces_per_cell); + storage_length * ReferenceCell::max_n_faces()); face_data_by_cells[my_q] .jacobian_gradients_non_inverse[0] .resize_fast(storage_length * - GeometryInfo::faces_per_cell); + ReferenceCell::max_n_faces()); } if (update_flags & update_quadrature_points) face_data_by_cells[my_q].quadrature_points.resize_fast( - cell_type.size() * GeometryInfo::faces_per_cell * + cell_type.size() * ReferenceCell::max_n_faces() * 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::faces_per_cell + + .data_index_offsets[cell * ReferenceCell::max_n_faces() + 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::faces_per_cell + face] + + [cell * ReferenceCell::max_n_faces() + 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::faces_per_cell * sizeof(GeometryType); + ReferenceCell::max_n_faces() * 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::faces_per_cell * + ReferenceCell::max_n_faces() * sizeof(GeometryType)); for (unsigned int j = 0; j < cell_data.size(); ++j) diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index a68bc1be24..4eecd5539a 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -2571,7 +2571,7 @@ MatrixFree::get_faces_by_cells_boundary_id( const unsigned int face_number) const { AssertIndexRange(cell_batch_index, n_cell_batches()); - AssertIndexRange(face_number, GeometryInfo::faces_per_cell); + AssertIndexRange(face_number, ReferenceCell::max_n_faces()); Assert(face_info.cell_and_face_boundary_id.size(0) >= n_cell_batches(), ExcNotInitialized()); std::array result; diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 11c72fd1b6..23affe7c90 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -2062,7 +2062,7 @@ MatrixFree::initialize_indices( numbers::invalid_unsigned_int); face_info.cell_and_face_boundary_id.reinit( TableIndices<3>(task_info.cell_partition_data.back(), - GeometryInfo::faces_per_cell, + ReferenceCell::max_n_faces(), VectorizedArrayType::size()), true); face_info.cell_and_face_boundary_id.fill(numbers::invalid_boundary_id); diff --git a/source/grid/grid_in.cc b/source/grid/grid_in.cc index 01523b6140..909902b64e 100644 --- a/source/grid/grid_in.cc +++ b/source/grid/grid_in.cc @@ -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::faces_per_cell; std::map> 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() + 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 &cell = cells[face_id / max_faces_per_cell]; - const ReferenceCell cell_type = + const unsigned int local_face_n = + face_id % ReferenceCell::max_n_faces(); + const CellData &cell = + cells[face_id / ReferenceCell::max_n_faces()]; + 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); diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 02305f7aee..6eacddc97c 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -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(-1)) - level.neighbors[cell * GeometryInfo::faces_per_cell + + level.neighbors[cell * ReferenceCell::max_n_faces() + j] = {0, nei.col[i]}; // set face indices - cells_0.cells[cell * GeometryInfo::faces_per_cell + j] = + cells_0.cells[cell * ReferenceCell::max_n_faces() + j] = crs.col[i]; // set face orientation if needed if (orientation_needed) { level.face_orientations.set_combined_orientation( - cell * GeometryInfo::faces_per_cell + j, + cell * ReferenceCell::max_n_faces() + j, connectivity.entity_orientations(dim - 1) .get_combined_orientation(i)); } @@ -15976,7 +15974,7 @@ void Triangulation::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() + face); std::array vertex_order{ {ref_cell.standard_to_real_face_vertex(0, face, diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 7a2bf6f9c4..704ae97b04 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -2360,19 +2360,19 @@ CellAccessor::set_neighbor( if (pointer.state() == IteratorState::valid) { this->tria->levels[this->present_level] - ->neighbors[this->present_index * GeometryInfo::faces_per_cell + i] + ->neighbors[this->present_index * ReferenceCell::max_n_faces() + i] .first = pointer->present_level; this->tria->levels[this->present_level] - ->neighbors[this->present_index * GeometryInfo::faces_per_cell + i] + ->neighbors[this->present_index * ReferenceCell::max_n_faces() + i] .second = pointer->present_index; } else { this->tria->levels[this->present_level] - ->neighbors[this->present_index * GeometryInfo::faces_per_cell + i] + ->neighbors[this->present_index * ReferenceCell::max_n_faces() + i] .first = -1; this->tria->levels[this->present_level] - ->neighbors[this->present_index * GeometryInfo::faces_per_cell + i] + ->neighbors[this->present_index * ReferenceCell::max_n_faces() + i] .second = -1; } } -- 2.39.5