From 4e22303d2bef71a608a4ac69f7c1e695fad1bc16 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sun, 24 Mar 2024 16:34:27 -0400 Subject: [PATCH] PeriodicFacePair: Use unsigned char instead of bitset. In addition, use the standard (orientation, rotation, flip) order in a temporary translation step. --- .../incompatibilities/20240525DavidWells | 8 ++ include/deal.II/dofs/dof_tools.h | 105 ++---------------- include/deal.II/grid/grid_tools.h | 60 ++-------- include/deal.II/grid/tria.h | 4 +- .../deal.II/matrix_free/face_setup_internal.h | 25 ++--- source/distributed/tria.cc | 69 ++++++------ source/dofs/dof_tools_constraints.cc | 80 +++++-------- source/dofs/dof_tools_constraints.inst.in | 12 +- source/grid/grid_tools.cc | 8 +- source/grid/grid_tools_dof_handlers.cc | 23 ++-- source/grid/grid_tools_dof_handlers.inst.in | 4 +- source/grid/tria.cc | 35 +++--- source/grid/tria_accessor.cc | 10 +- source/multigrid/mg_constrained_dofs.cc | 4 +- tests/dofs/dof_tools_21_b.cc | 27 +++-- tests/dofs/dof_tools_21_b_x.cc | 15 +-- tests/dofs/dof_tools_21_b_x_q3.cc | 15 +-- tests/dofs/dof_tools_21_b_y.cc | 15 +-- tests/dofs/dof_tools_21_c.cc | 27 +++-- tests/grid/grid_tools_05.cc | 12 +- tests/grid/grid_tools_06.cc | 12 +- tests/grid/periodicity_1d.cc | 7 +- tests/mpi/periodicity_04.cc | 6 +- 23 files changed, 193 insertions(+), 390 deletions(-) create mode 100644 doc/news/changes/incompatibilities/20240525DavidWells diff --git a/doc/news/changes/incompatibilities/20240525DavidWells b/doc/news/changes/incompatibilities/20240525DavidWells new file mode 100644 index 0000000000..bbbbb75e43 --- /dev/null +++ b/doc/news/changes/incompatibilities/20240525DavidWells @@ -0,0 +1,8 @@ +Changed: GridTools::orthogonal_equality() now returns a +std::optional instead of a +std::optional>. As a consequence, +PeriodicFacePair now also stores an unsigned char instead of +a std::bitset<3> to represent the relative orientation of the first face to the +second face. +
+(David Wells, 2024/05/25) diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index 884bf558ce..bb380b6440 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -917,97 +917,9 @@ namespace DoFTools * of components of the finite element. This can be used to enforce * periodicity in only one variable in a system of equations. * - * @p face_orientation, @p face_flip and @p face_rotation describe an - * orientation that should be applied to @p face_1 prior to matching and - * constraining DoFs. This has nothing to do with the actual orientation of - * the given faces in their respective cells (which for boundary faces is - * always the default) but instead how you want to see periodicity to be - * enforced. For example, by using these flags, you can enforce a condition - * of the kind $u(0,y)=u(1,1-y)$ (i.e., a Moebius band) or in 3d a twisted - * torus. More precisely, these flags match local face DoF indices in the - * following manner: - * - * In 2d: face_orientation must always be true, - * face_rotation is always false, and face_flip has the - * meaning of line_flip; this implies e.g. for Q1: - * - * @code - * - * face_orientation = true, face_flip = false, face_rotation = false: - * - * face1: face2: - * - * 1 1 - * | <--> | - * 0 0 - * - * Resulting constraints: 0 <-> 0, 1 <-> 1 - * - * (Numbers denote local face DoF indices.) - * - * - * face_orientation = true, face_flip = true, face_rotation = false: - * - * face1: face2: - * - * 0 1 - * | <--> | - * 1 0 - * - * Resulting constraints: 1 <-> 0, 0 <-> 1 - * @endcode - * - * And similarly for the case of Q1 in 3d: - * - * @code - * - * face_orientation = true, face_flip = false, face_rotation = false: - * - * face1: face2: - * - * 2 - 3 2 - 3 - * | | <--> | | - * 0 - 1 0 - 1 - * - * Resulting constraints: 0 <-> 0, 1 <-> 1, 2 <-> 2, 3 <-> 3 - * - * (Numbers denote local face DoF indices.) - * - * - * face_orientation = false, face_flip = false, face_rotation = false: - * - * face1: face2: - * - * 1 - 3 2 - 3 - * | | <--> | | - * 0 - 2 0 - 1 - * - * Resulting constraints: 0 <-> 0, 2 <-> 1, 1 <-> 2, 3 <-> 3 - * - * - * face_orientation = true, face_flip = true, face_rotation = false: - * - * face1: face2: - * - * 1 - 0 2 - 3 - * | | <--> | | - * 3 - 2 0 - 1 - * - * Resulting constraints: 3 <-> 0, 2 <-> 1, 1 <-> 2, 0 <-> 3 - * - * - * face_orientation = true, face_flip = false, face_rotation = true - * - * face1: face2: - * - * 0 - 2 2 - 3 - * | | <--> | | - * 1 - 3 0 - 1 - * - * Resulting constraints: 1 <-> 0, 3 <-> 1, 0 <-> 2, 2 <-> 3 - * - * and any combination of that... - * @endcode + * Here, @p combined_orientation is the relative orientation of @p face_1 with + * respect to @p face_2. This is typically computed by + * GridTools::orthogonal_equality(). * * Optionally a matrix @p matrix along with a std::vector @p * first_vector_components can be specified that describes how DoFs on @p @@ -1045,10 +957,9 @@ namespace DoFTools const FaceIterator &face_1, const std_cxx20::type_identity_t &face_2, AffineConstraints &constraints, - const ComponentMask &component_mask = {}, - const bool face_orientation = true, - const bool face_flip = false, - const bool face_rotation = false, + const ComponentMask &component_mask = {}, + const unsigned char combined_orientation = + ReferenceCell::default_combined_face_orientation(), const FullMatrix &matrix = FullMatrix(), const std::vector &first_vector_components = std::vector(), @@ -1198,9 +1109,7 @@ namespace DoFTools const FullMatrix &transformation, AffineConstraints &affine_constraints, const ComponentMask &component_mask, - const bool face_orientation, - const bool face_flip, - const bool face_rotation, + const unsigned char combined_orientation, const number periodicity_factor, const unsigned int level = numbers::invalid_unsigned_int); } // namespace internal diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 04b19802a3..5fc2e5dfd0 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -58,7 +58,6 @@ # include #endif -#include #include #include @@ -2336,9 +2335,9 @@ namespace GridTools /** * The relative orientation of the first face with respect to the second * face as described in orthogonal_equality() and - * DoFTools::make_periodicity_constraints() (and stored as a bitset). + * DoFTools::make_periodicity_constraints(). */ - std::bitset<3> orientation; + unsigned char orientation; /** * A @p dim $\times$ @p dim rotation matrix that describes how vector @@ -2377,59 +2376,14 @@ namespace GridTools * identity matrix. * * If the matching was successful, the _relative_ orientation of @p face1 with - * respect to @p face2 is returned a std::optional> object - * orientation in which - * @code - * orientation.value()[0] = face_orientation - * orientation.value()[1] = face_flip - * orientation.value()[2] = face_rotation - * @endcode - * - * In 2d face_orientation is always true, - * face_rotation is always false, and face_flip has the - * meaning of line_flip. More precisely in 3d: - * - * face_orientation: true if @p face1 and @p face2 have - * the same orientation. Otherwise, the vertex indices of @p face1 match the - * vertex indices of @p face2 in the following manner: - * - * @code - * face1: face2: - * - * 1 - 3 2 - 3 - * | | <--> | | - * 0 - 2 0 - 1 - * @endcode - * - * face_flip: true if the matched vertices are rotated by - * 180 degrees: - * - * @code - * face1: face2: - * - * 1 - 0 2 - 3 - * | | <--> | | - * 3 - 2 0 - 1 - * @endcode - * - * face_rotation: true if the matched vertices are rotated - * by 90 degrees counterclockwise: - * - * @code - * face1: face2: - * - * 0 - 2 2 - 3 - * | | <--> | | - * 1 - 3 0 - 1 - * @endcode - * - * and any combination of that... More information on the topic can be found - * in the + * respect to @p face2 is returned a std::optional, in which + * the stored value is the same orientation bit format used elsewhere in the + * library. More information on that topic can be found in the * @ref GlossFaceOrientation "glossary" * article. */ template - std::optional> + std::optional orthogonal_equality( const FaceIterator &face1, const FaceIterator &face2, @@ -2451,7 +2405,7 @@ namespace GridTools * with faces belonging to the second boundary with the help of * orthogonal_equality(). * - * The bitset that is returned inside of PeriodicFacePair encodes the + * The unsigned char that is returned inside of PeriodicFacePair encodes the * _relative_ orientation of the first face with respect to the second face, * see the documentation of orthogonal_equality() for further details. * diff --git a/include/deal.II/grid/tria.h b/include/deal.II/grid/tria.h index 3118a3ed23..58e08424b2 100644 --- a/include/deal.II/grid/tria.h +++ b/include/deal.II/grid/tria.h @@ -3649,7 +3649,7 @@ public: */ const std::map< std::pair, - std::pair, std::bitset<3>>> & + std::pair, unsigned char>> & get_periodic_face_map() const; /** @@ -4089,7 +4089,7 @@ private: * face pairs. */ std::map, - std::pair, std::bitset<3>>> + std::pair, unsigned char>> periodic_face_map; /** diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h index 3a03f5bec7..453d1c7571 100644 --- a/include/deal.II/matrix_free/face_setup_internal.h +++ b/include/deal.II/matrix_free/face_setup_internal.h @@ -1004,21 +1004,16 @@ namespace internal // special treatment of periodic boundaries if (dim == 3 && cell->has_periodic_neighbor(face_no)) { - const unsigned int exterior_face_orientation = - !cell->get_triangulation() - .get_periodic_face_map() - .at({cell, face_no}) - .second[0] + - 2 * cell->get_triangulation() - .get_periodic_face_map() - .at({cell, face_no}) - .second[1] + - 4 * cell->get_triangulation() - .get_periodic_face_map() - .at({cell, face_no}) - .second[2]; - - info.face_orientation = exterior_face_orientation; + unsigned char exterior_face_orientation = cell->get_triangulation() + .get_periodic_face_map() + .at({cell, face_no}) + .second; + const auto [orientation, rotation, flip] = + ::dealii::internal::split_face_orientation( + exterior_face_orientation); + + info.face_orientation = + (orientation ? 0u : 1u) + 2 * flip + 4 * rotation; return info; } diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc index 753c7f7ebc..651349c8e2 100644 --- a/source/distributed/tria.cc +++ b/source/distributed/tria.cc @@ -2590,16 +2590,18 @@ namespace parallel typename Triangulation::cell_iterator; typename std::map, std::pair, - std::bitset<3>>>::const_iterator it; + unsigned char>>::const_iterator it; for (it = tria.get_periodic_face_map().begin(); it != tria.get_periodic_face_map().end(); ++it) { - const cell_iterator &cell_1 = it->first.first; - const unsigned int face_no_1 = it->first.second; - const cell_iterator &cell_2 = it->second.first.first; - const unsigned int face_no_2 = it->second.first.second; - const std::bitset<3> face_orientation = it->second.second; + const cell_iterator &cell_1 = it->first.first; + const unsigned int face_no_1 = it->first.second; + const cell_iterator &cell_2 = it->second.first.first; + const unsigned int face_no_2 = it->second.first.second; + const unsigned char combined_orientation = it->second.second; + const auto [orientation, rotation, flip] = + ::dealii::internal::split_face_orientation(combined_orientation); if (cell_1->level() == cell_2->level()) { @@ -2611,10 +2613,7 @@ namespace parallel // cell[0] into account const unsigned int vface0 = GeometryInfo::standard_to_real_face_vertex( - v, - face_orientation[0], - face_orientation[1], - face_orientation[2]); + v, orientation, flip, rotation); const unsigned int vi0 = topological_vertex_numbering[cell_1->face(face_no_1) ->vertex_index(vface0)]; @@ -3706,11 +3705,13 @@ namespace parallel for (unsigned int repetition = 0; repetition < dim; ++repetition) for (const auto &it : this->get_periodic_face_map()) { - const cell_iterator &cell_1 = it.first.first; - const unsigned int face_no_1 = it.first.second; - const cell_iterator &cell_2 = it.second.first.first; - const unsigned int face_no_2 = it.second.first.second; - const std::bitset<3> &face_orientation = it.second.second; + const cell_iterator &cell_1 = it.first.first; + const unsigned int face_no_1 = it.first.second; + const cell_iterator &cell_2 = it.second.first.first; + const unsigned int face_no_2 = it.second.first.second; + const unsigned char combined_orientation = it.second.second; + const auto [orientation, rotation, flip] = + ::dealii::internal::split_face_orientation(combined_orientation); if (cell_1->level() == level && cell_2->level() == level) { @@ -3722,10 +3723,7 @@ namespace parallel // account const unsigned int vface0 = GeometryInfo::standard_to_real_face_vertex( - v, - face_orientation[0], - face_orientation[1], - face_orientation[2]); + v, orientation, flip, rotation); if (marked_vertices[cell_1->face(face_no_1)->vertex_index( vface0)] || marked_vertices[cell_2->face(face_no_2)->vertex_index( @@ -3794,16 +3792,22 @@ namespace parallel const unsigned int tree_right = coarse_cell_to_p4est_tree_permutation[second_cell->index()]; - // p4est wants to know which corner the first corner on - // the face with the lower id is mapped to on the face with - // with the higher id. For d==2 there are only two possibilities - // that are determined by it->orientation[0] (the face_orientation - // flag). For d==3 we have to use GridTools::OrientationLookupTable. - // The result is given below. + // p4est wants to know which corner the first corner on the face with + // the lower id is mapped to on the face with with the higher id. For + // d==2 there are only two possibilities: i.e., face_pair.orientation + // must be 0 or 1. For d==3 we have to use a lookup table. The result + // is given below. unsigned int p4est_orientation = 0; if (dim == 2) - p4est_orientation = face_pair.orientation[0] == true ? 0u : 1u; + { + AssertIndexRange(face_pair.orientation, 2); + p4est_orientation = + face_pair.orientation == + ReferenceCell::default_combined_face_orientation() ? + 0u : + 1u; + } else { const unsigned int face_idx_list[] = {face_left, face_right}; @@ -3847,26 +3851,27 @@ namespace parallel Assert(first_dealii_idx_on_face != numbers::invalid_unsigned_int, ExcInternalError()); // Now map dealii_idx_on_face according to the orientation + constexpr unsigned int left_to_right[8][4] = {{0, 2, 1, 3}, {0, 1, 2, 3}, - {3, 1, 2, 0}, - {3, 2, 1, 0}, {2, 3, 0, 1}, {1, 3, 0, 2}, + {3, 1, 2, 0}, + {3, 2, 1, 0}, {1, 0, 3, 2}, {2, 0, 3, 1}}; constexpr unsigned int right_to_left[8][4] = {{0, 2, 1, 3}, {0, 1, 2, 3}, - {3, 1, 2, 0}, - {3, 2, 1, 0}, {2, 3, 0, 1}, {2, 0, 3, 1}, + {3, 1, 2, 0}, + {3, 2, 1, 0}, {1, 0, 3, 2}, {1, 3, 0, 2}}; const unsigned int second_dealii_idx_on_face = - lower_idx == 0 ? left_to_right[face_pair.orientation.to_ulong()] + lower_idx == 0 ? left_to_right[face_pair.orientation] [first_dealii_idx_on_face] : - right_to_left[face_pair.orientation.to_ulong()] + right_to_left[face_pair.orientation] [first_dealii_idx_on_face]; const unsigned int second_dealii_idx_on_cell = GeometryInfo::face_to_cell_vertices( diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index f15482d330..56f32a40bc 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -3106,20 +3106,13 @@ namespace DoFTools const FullMatrix &transformation, AffineConstraints &affine_constraints, const ComponentMask &component_mask, - const bool face_orientation, - const bool face_flip, - const bool face_rotation, + const unsigned char combined_orientation, const number periodicity_factor, const unsigned int level) { static const int dim = FaceIterator::AccessorType::dimension; static const int spacedim = FaceIterator::AccessorType::space_dimension; - const unsigned char combined_face_orientation = - ::dealii::internal::combined_face_orientation(face_orientation, - face_rotation, - face_flip); - const bool use_mg = (level != numbers::invalid_unsigned_int); // If we don't use multigrid, we should be in the case where face_1 is @@ -3174,9 +3167,7 @@ namespace DoFTools child_transformation, affine_constraints, component_mask, - face_orientation, - face_flip, - face_rotation, + combined_orientation, periodicity_factor); } return; @@ -3362,7 +3353,7 @@ namespace DoFTools // given orientation: const unsigned int j = cell_to_rotated_face_index[fe.face_to_cell_index( - jj, 0, combined_face_orientation)]; + jj, 0, combined_orientation)]; if (std::abs(transformation(i, jj)) > eps) constraint_entries.emplace_back(dofs_1[j], @@ -3385,7 +3376,7 @@ namespace DoFTools // orientation: const unsigned int j = cell_to_rotated_face_index[fe.face_to_cell_index( - target, 0, combined_face_orientation)]; + target, 0, combined_orientation)]; auto dof_left = dofs_1[j]; auto dof_right = dofs_2[i]; @@ -3601,9 +3592,7 @@ namespace DoFTools const std_cxx20::type_identity_t &face_2, AffineConstraints &affine_constraints, const ComponentMask &component_mask, - const bool face_orientation, - const bool face_flip, - const bool face_rotation, + const unsigned char combined_orientation, const FullMatrix &matrix, const std::vector &first_vector_components, const number periodicity_factor) @@ -3611,15 +3600,17 @@ namespace DoFTools static const int dim = FaceIterator::AccessorType::dimension; static const int spacedim = FaceIterator::AccessorType::space_dimension; - Assert((dim != 1) || (face_orientation == true && face_flip == false && - face_rotation == false), - ExcMessage("The supplied orientation " - "(face_orientation, face_flip, face_rotation) " +#ifdef DEBUG + const auto [orientation, rotation, flip] = + ::dealii::internal::split_face_orientation(combined_orientation); + + Assert((dim != 1) || + (orientation == true && flip == false && rotation == false), + ExcMessage("The supplied orientation (orientation, rotation, flip) " "is invalid for 1d")); - Assert((dim != 2) || (face_flip == false && face_rotation == false), - ExcMessage("The supplied orientation " - "(face_orientation, face_flip, face_rotation) " + Assert((dim != 2) || (flip == false && rotation == false), + ExcMessage("The supplied orientation (orientation, rotation, flip) " "is invalid for 2d")); Assert(face_1 != face_2, @@ -3638,12 +3629,6 @@ namespace DoFTools ExcMessage("first_vector_components is nonempty, so matrix must " "be a rotation matrix exactly of size spacedim")); - const unsigned char combined_orientation = - ::dealii::internal::combined_face_orientation(face_orientation, - face_rotation, - face_flip); - -#ifdef DEBUG if (!face_1->has_children()) { // TODO: the implementation makes the assumption that all faces have the @@ -3724,9 +3709,7 @@ namespace DoFTools face_2->child(j), affine_constraints, component_mask, - face_orientation, - face_flip, - face_rotation, + combined_orientation, matrix, first_vector_components, periodicity_factor); @@ -3769,9 +3752,7 @@ namespace DoFTools transformation, affine_constraints, component_mask, - face_orientation, - face_flip, - face_rotation, + combined_orientation, periodicity_factor); } else @@ -3784,9 +3765,7 @@ namespace DoFTools inverse, affine_constraints, component_mask, - face_orientation, - face_flip, - face_rotation, + combined_orientation, periodicity_factor); } } @@ -3800,17 +3779,16 @@ namespace DoFTools // the rotation when constraining face_2 to face_1. Therefore // face_flip has to be toggled if face_rotation is true: In case of // inverted orientation, nothing has to be done. - internal::set_periodicity_constraints(face_1, - face_2, - transformation, - affine_constraints, - component_mask, - face_orientation, - face_orientation ? - face_rotation ^ face_flip : - face_flip, - face_rotation, - periodicity_factor); + + internal::set_periodicity_constraints( + face_1, + face_2, + transformation, + affine_constraints, + component_mask, + ::dealii::internal::combined_face_orientation( + orientation, rotation, orientation ? rotation ^ flip : flip), + periodicity_factor); } } } @@ -3845,9 +3823,7 @@ namespace DoFTools face_2, constraints, component_mask, - pair.orientation[0], - pair.orientation[1], - pair.orientation[2], + pair.orientation, pair.matrix, first_vector_components, periodicity_factor); diff --git a/source/dofs/dof_tools_constraints.inst.in b/source/dofs/dof_tools_constraints.inst.in index c2746d5e8e..c13b320141 100644 --- a/source/dofs/dof_tools_constraints.inst.in +++ b/source/dofs/dof_tools_constraints.inst.in @@ -23,9 +23,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; deal_II_space_dimension>::face_iterator &, AffineConstraints &, const ComponentMask &, - bool, - bool, - bool, + const unsigned char, const FullMatrix &, const std::vector &, const S); @@ -38,9 +36,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; const FullMatrix &transformation, AffineConstraints &affine_constraints, const ComponentMask &component_mask, - const bool face_orientation, - const bool face_flip, - const bool face_rotation, + const unsigned char combined_orientation, const S periodicity_factor, const unsigned int level); @@ -52,9 +48,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; const FullMatrix &transformation, AffineConstraints &affine_constraints, const ComponentMask &component_mask, - const bool face_orientation, - const bool face_flip, - const bool face_rotation, + const unsigned char combined_orientation, const S periodicity_factor, const unsigned int level); #endif diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 5b84682905..edb9ee540a 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -4826,16 +4826,16 @@ namespace GridTools const auto face_a = pair.first.first->face(pair.first.second); const auto face_b = pair.second.first.first->face(pair.second.first.second); - const auto mask = pair.second.second; + const unsigned char combined_orientation = pair.second.second; AssertDimension(face_a->n_vertices(), face_b->n_vertices()); // loop over all vertices on face for (unsigned int i = 0; i < face_a->n_vertices(); ++i) { - const bool face_orientation = mask[0]; - const bool face_flip = mask[1]; - const bool face_rotation = mask[2]; + const auto [face_orientation, face_rotation, face_flip] = + ::dealii::internal::split_face_orientation( + combined_orientation); // find the right local vertex index for the second face unsigned int j = 0; diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index 8d9fb293f7..f7cdc8f9f2 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -2163,7 +2163,7 @@ namespace GridTools const CellIterator cell2 = it2->first; const unsigned int face_idx1 = it1->second; const unsigned int face_idx2 = it2->second; - if (const std::optional> orientation = + if (const std::optional orientation = GridTools::orthogonal_equality(cell1->face(face_idx1), cell2->face(face_idx2), direction, @@ -2419,7 +2419,7 @@ namespace GridTools template - std::optional> + std::optional orthogonal_equality( const FaceIterator &face1, const FaceIterator &face2, @@ -2472,7 +2472,7 @@ namespace GridTools face1_vertices.cbegin() + face1->n_vertices()), make_array_view(face2_vertices.cbegin(), face2_vertices.cbegin() + face2->n_vertices())); - std::bitset<3> orientation; + unsigned char orientation = std::numeric_limits::max(); if (dim == 1) { // In 1D things are always well-oriented @@ -2489,20 +2489,11 @@ namespace GridTools else { Assert(dim == 3, ExcInternalError()); - // There are two differences between the orientation implementation - // used in the periodicity code and that used in ReferenceCell: - // - // 1. The bitset is unpacked as (orientation, flip, rotation) - // instead of the standard (orientation, rotation, flip). - // - // 2. The 90 degree rotations are always clockwise, so the third and - // seventh (in the combined orientation) are switched. - // - // Both translations are encoded in this table. This matches - // OrientationLookupTable<3> which was present in previous revisions - // of this file. + // Unlike the standard orientation, here the 90 degree rotations are + // always clockwise, so the third and seventh (in the combined + // orientation) are switched. constexpr std::array translation{ - {0, 1, 4, 7, 2, 3, 6, 5}}; + {0, 1, 2, 7, 4, 5, 6, 3}}; AssertIndexRange(combined_orientation, translation.size()); orientation = translation[std::min(combined_orientation, 7u)]; diff --git a/source/grid/grid_tools_dof_handlers.inst.in b/source/grid/grid_tools_dof_handlers.inst.in index c53ff46a0f..e7599cc7e5 100644 --- a/source/grid/grid_tools_dof_handlers.inst.in +++ b/source/grid/grid_tools_dof_handlers.inst.in @@ -249,7 +249,7 @@ for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLER; namespace GridTools \{ - template std::optional> + template std::optional orthogonal_equality( const X::active_face_iterator &, const X::active_face_iterator &, @@ -257,7 +257,7 @@ for (X : SEQUENTIAL_TRIANGULATION_AND_DOFHANDLER; const Tensor<1, deal_II_space_dimension> &, const FullMatrix &); - template std::optional> + template std::optional orthogonal_equality( const X::face_iterator &, const X::face_iterator &, diff --git a/source/grid/tria.cc b/source/grid/tria.cc index 8d68212f98..bfd5252b33 100644 --- a/source/grid/tria.cc +++ b/source/grid/tria.cc @@ -1643,21 +1643,20 @@ namespace const typename Triangulation::cell_iterator &cell_2, unsigned int n_face_1, unsigned int n_face_2, - const std::bitset<3> &orientation, + const unsigned char orientation, typename std::map< std::pair::cell_iterator, unsigned int>, std::pair::cell_iterator, unsigned int>, - std::bitset<3>>> &periodic_face_map) + unsigned char>> &periodic_face_map) { using FaceIterator = typename Triangulation::face_iterator; const FaceIterator face_1 = cell_1->face(n_face_1); const FaceIterator face_2 = cell_2->face(n_face_2); - const bool face_orientation = orientation[0]; - const bool face_flip = orientation[1]; - const bool face_rotation = orientation[2]; + const auto [face_orientation, face_rotation, face_flip] = + internal::split_face_orientation(orientation); Assert((dim != 1) || (face_orientation == true && face_flip == false && face_rotation == false), @@ -1686,12 +1685,12 @@ namespace using CellFace = std::pair::cell_iterator, unsigned int>; - const CellFace cell_face_1(cell_1, n_face_1); - const CellFace cell_face_2(cell_2, n_face_2); - const std::pair> cell_face_orientation_2( + const CellFace cell_face_1(cell_1, n_face_1); + const CellFace cell_face_2(cell_2, n_face_2); + const std::pair cell_face_orientation_2( cell_face_2, orientation); - const std::pair>> + const std::pair> periodic_faces(cell_face_1, cell_face_orientation_2); // Only one periodic neighbor is allowed @@ -15790,7 +15789,7 @@ const typename std::map< std::pair::cell_iterator, unsigned int>, std::pair::cell_iterator, unsigned int>, - std::bitset<3>>> + unsigned char>> &Triangulation::get_periodic_face_map() const { return periodic_face_map; @@ -16027,15 +16026,13 @@ void Triangulation::update_periodic_face_map() periodic_face_map); // for the other way, we need to invert the orientation - std::bitset<3> inverted_orientation; + unsigned char inverted_orientation; { - bool orientation, flip, rotation; - orientation = it->orientation[0]; - rotation = it->orientation[2]; - flip = orientation ? rotation ^ it->orientation[1] : it->orientation[1]; - inverted_orientation[0] = orientation; - inverted_orientation[1] = flip; - inverted_orientation[2] = rotation; + auto [orientation, rotation, flip] = + internal::split_face_orientation(it->orientation); + flip = orientation ? rotation ^ flip : flip; + inverted_orientation = + internal::combined_face_orientation(orientation, rotation, flip); } update_periodic_face_map_recursively(it->cell[1], it->cell[0], @@ -16048,7 +16045,7 @@ void Triangulation::update_periodic_face_map() // check consistency typename std::map, std::pair, - std::bitset<3>>>::const_iterator it_test; + unsigned char>>::const_iterator it_test; for (it_test = periodic_face_map.begin(); it_test != periodic_face_map.end(); ++it_test) { diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index a6522524ce..000a1c2158 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -2801,13 +2801,17 @@ CellAccessor::periodic_neighbor_child_on_subface( * number of children as i_subface. */ AssertIndexRange(i_subface, nb_parent_face_it->n_children()); + + const auto [orientation, rotation, flip] = + internal::split_face_orientation(my_face_pair->second.second); + unsigned int sub_neighbor_num = GeometryInfo::child_cell_on_face(parent_nb_it->refinement_case(), nb_face_num, i_subface, - my_face_pair->second.second[0], - my_face_pair->second.second[1], - my_face_pair->second.second[2], + orientation, + flip, + rotation, nb_parent_face_it->refinement_case()); return parent_nb_it->child(sub_neighbor_num); } diff --git a/source/multigrid/mg_constrained_dofs.cc b/source/multigrid/mg_constrained_dofs.cc index d9fccb66f7..c1effb4ff3 100644 --- a/source/multigrid/mg_constrained_dofs.cc +++ b/source/multigrid/mg_constrained_dofs.cc @@ -88,9 +88,7 @@ MGConstrainedDoFs::initialize( transformation, level_constraints[first_cell.first->level()], component_mask, - second_cell.second[0], - second_cell.second[1], - second_cell.second[2], + second_cell.second, periodicity_factor, first_cell.first->level()); } diff --git a/tests/dofs/dof_tools_21_b.cc b/tests/dofs/dof_tools_21_b.cc index cfdafcd533..03d1928332 100644 --- a/tests/dofs/dof_tools_21_b.cc +++ b/tests/dofs/dof_tools_21_b.cc @@ -26,6 +26,7 @@ #include +#include #include #include @@ -34,11 +35,11 @@ // // Test // DoFTools:: -// make_periodicity_constraints (const FaceIterator &, -// const FaceIterator &, -// dealii::AffineConstraints &, -// const std::vector &, -// bool, bool, bool) +// make_periodicity_constraints(const FaceIterator &, +// const FaceIterator &, +// dealii::AffineConstraints &, +// const std::vector &, +// unsigned char) // for correct behavior on non standard oriented meshes. // @@ -320,10 +321,13 @@ print_matching(DoFHandler &dof_handler, face_1, face_2, dim == 2 ? 1 : 2, Tensor<1, spacedim>()); AssertThrow(orientation, ExcInternalError()); const auto o = *orientation; - deallog << "Orientation: " << o[0] << o[1] << o[2] << std::endl; + // Preserve old output by printing a bitset and also using the old + // (orientation, flip, rotation) format + std::bitset<3> bo(o); + deallog << "Orientation: " << bo[0] << bo[2] << bo[1] << std::endl; DoFTools::make_periodicity_constraints( - face_1, face_2, constraint_matrix, velocity_mask, o[0], o[1], o[2]); + face_1, face_2, constraint_matrix, velocity_mask, o); deallog << "Matching:" << std::endl; constraint_matrix.print(deallog.get_file_stream()); constraint_matrix.close(); @@ -332,13 +336,8 @@ print_matching(DoFHandler &dof_handler, face_2, face_1, dim == 2 ? 1 : 2, Tensor<1, spacedim>()); AssertThrow(reverse_orientation, ExcInternalError()); const auto ro = *reverse_orientation; - DoFTools::make_periodicity_constraints(face_2, - face_1, - constraint_matrix_reverse, - velocity_mask, - ro[0], - ro[1], - ro[2]); + DoFTools::make_periodicity_constraints( + face_2, face_1, constraint_matrix_reverse, velocity_mask, ro); deallog << "Reverse Matching:" << std::endl; constraint_matrix_reverse.print(deallog.get_file_stream()); constraint_matrix_reverse.close(); diff --git a/tests/dofs/dof_tools_21_b_x.cc b/tests/dofs/dof_tools_21_b_x.cc index 74b04cae30..228a524c86 100644 --- a/tests/dofs/dof_tools_21_b_x.cc +++ b/tests/dofs/dof_tools_21_b_x.cc @@ -195,19 +195,8 @@ print_matching(DoFHandler &dof_handler) deallog << dofs_2[i] << " is located at " << support_points[dofs_2[i]] << std::endl; - - std::bitset<3> orientation; - orientation[0] = 0; - orientation[1] = 0; - orientation[2] = 0; - - DoFTools::make_periodicity_constraints(face_1, - face_2, - constraint_matrix, - ComponentMask(), - orientation[0], - orientation[1], - orientation[2]); + DoFTools::make_periodicity_constraints( + face_1, face_2, constraint_matrix, ComponentMask(), 0u); constraint_matrix.print(deallog.get_file_stream()); constraint_matrix.close(); deallog << "Matching:" << std::endl; diff --git a/tests/dofs/dof_tools_21_b_x_q3.cc b/tests/dofs/dof_tools_21_b_x_q3.cc index 96f212b80c..245a2bcc28 100644 --- a/tests/dofs/dof_tools_21_b_x_q3.cc +++ b/tests/dofs/dof_tools_21_b_x_q3.cc @@ -229,19 +229,8 @@ print_matching(DoFHandler &dof_handler) deallog << dofs_2[i] << " is located at " << support_points[dofs_2[i]] << std::endl; - - std::bitset<3> orientation; - orientation[0] = 0; - orientation[1] = 0; - orientation[2] = 0; - - DoFTools::make_periodicity_constraints(face_1, - face_2, - constraint_matrix, - ComponentMask(), - orientation[0], - orientation[1], - orientation[2]); + DoFTools::make_periodicity_constraints( + face_1, face_2, constraint_matrix, ComponentMask(), 0u); constraint_matrix.print(deallog.get_file_stream()); constraint_matrix.close(); deallog << "Matching:" << std::endl; diff --git a/tests/dofs/dof_tools_21_b_y.cc b/tests/dofs/dof_tools_21_b_y.cc index 485f24f1d3..1763051683 100644 --- a/tests/dofs/dof_tools_21_b_y.cc +++ b/tests/dofs/dof_tools_21_b_y.cc @@ -195,19 +195,8 @@ print_matching(DoFHandler &dof_handler) deallog << dofs_2[i] << " is located at " << support_points[dofs_2[i]] << std::endl; - - std::bitset<3> orientation; - orientation[0] = 0; - orientation[1] = 0; - orientation[2] = 0; - - DoFTools::make_periodicity_constraints(face_1, - face_2, - constraint_matrix, - ComponentMask(), - orientation[0], - orientation[1], - orientation[2]); + DoFTools::make_periodicity_constraints( + face_1, face_2, constraint_matrix, ComponentMask(), 0u); constraint_matrix.print(deallog.get_file_stream()); constraint_matrix.close(); deallog << "Matching:" << std::endl; diff --git a/tests/dofs/dof_tools_21_c.cc b/tests/dofs/dof_tools_21_c.cc index da89b9f026..3cee6d118e 100644 --- a/tests/dofs/dof_tools_21_c.cc +++ b/tests/dofs/dof_tools_21_c.cc @@ -26,6 +26,7 @@ #include +#include #include #include @@ -34,11 +35,11 @@ // // Test // DoFTools:: -// make_periodicity_constraints (const FaceIterator &, -// const FaceIterator &, -// dealii::AffineConstraints &, -// const std::vector &, -// bool, bool, bool) +// make_periodicity_constraints(const FaceIterator &, +// const FaceIterator &, +// dealii::AffineConstraints &, +// const std::vector &, +// unsigned char) // // for correct behavior with hanging nodes. This is done by additionally // refining the second cube once. Test that constraining face_1 -> face_2 @@ -326,10 +327,13 @@ print_matching(DoFHandler &dof_handler, face_1, face_2, dim == 2 ? 1 : 2, Tensor<1, spacedim>()); AssertThrow(orientation, ExcInternalError()); const auto o = *orientation; - deallog << "Orientation: " << o[0] << o[1] << o[2] << std::endl; + // Preserve old output by printing a bitset and also using the old + // (orientation, flip, rotation) format + std::bitset<3> bo(o); + deallog << "Orientation: " << bo[0] << bo[2] << bo[1] << std::endl; DoFTools::make_periodicity_constraints( - face_1, face_2, constraint_matrix, velocity_mask, o[0], o[1], o[2]); + face_1, face_2, constraint_matrix, velocity_mask, o); deallog << "Matching:" << std::endl; constraint_matrix.print(deallog.get_file_stream()); constraint_matrix.close(); @@ -338,13 +342,8 @@ print_matching(DoFHandler &dof_handler, face_2, face_1, dim == 2 ? 1 : 2, Tensor<1, spacedim>()); AssertThrow(reverse_orientation, ExcInternalError()); const auto ro = *reverse_orientation; - DoFTools::make_periodicity_constraints(face_2, - face_1, - constraint_matrix_reverse, - velocity_mask, - ro[0], - ro[1], - ro[2]); + DoFTools::make_periodicity_constraints( + face_2, face_1, constraint_matrix_reverse, velocity_mask, ro); deallog << "Reverse Matching:" << std::endl; constraint_matrix_reverse.print(deallog.get_file_stream()); constraint_matrix_reverse.close(); diff --git a/tests/grid/grid_tools_05.cc b/tests/grid/grid_tools_05.cc index 346b7e8a32..4e40aec722 100644 --- a/tests/grid/grid_tools_05.cc +++ b/tests/grid/grid_tools_05.cc @@ -154,9 +154,9 @@ generate_grid(Triangulation<3> &triangulation) */ template void -print_match(const FaceIterator &face_1, - const FaceIterator &face_2, - const std::bitset<3> &orientation) +print_match(const FaceIterator &face_1, + const FaceIterator &face_2, + const unsigned char combined_orientation) { static const int dim = FaceIterator::AccessorType::dimension; @@ -170,8 +170,10 @@ print_match(const FaceIterator &face_1, deallog << " :: " << face_2->vertex(j); deallog << std::endl; - deallog << "orientation: " << orientation[0] << " flip: " << orientation[1] - << " rotation: " << orientation[2] << std::endl + const auto [orientation, rotation, flip] = + internal::split_face_orientation(combined_orientation); + deallog << "orientation: " << orientation << " flip: " << flip + << " rotation: " << rotation << std::endl << std::endl; } diff --git a/tests/grid/grid_tools_06.cc b/tests/grid/grid_tools_06.cc index 2b222c990d..8ecfcfe76b 100644 --- a/tests/grid/grid_tools_06.cc +++ b/tests/grid/grid_tools_06.cc @@ -167,9 +167,9 @@ generate_grid(Triangulation<3> &triangulation, int orientation) */ template void -print_match(const FaceIterator &face_1, - const FaceIterator &face_2, - const std::bitset<3> &orientation) +print_match(const FaceIterator &face_1, + const FaceIterator &face_2, + const unsigned char combined_orientation) { static const int dim = FaceIterator::AccessorType::dimension; @@ -183,8 +183,10 @@ print_match(const FaceIterator &face_1, deallog << " :: " << face_2->vertex(j); deallog << std::endl; - deallog << "orientation: " << orientation[0] << " flip: " << orientation[1] - << " rotation: " << orientation[2] << std::endl + const auto [orientation, rotation, flip] = + internal::split_face_orientation(combined_orientation); + deallog << "orientation: " << orientation << " flip: " << flip + << " rotation: " << rotation << std::endl << std::endl; } diff --git a/tests/grid/periodicity_1d.cc b/tests/grid/periodicity_1d.cc index 006547a5fc..f21b2ae811 100644 --- a/tests/grid/periodicity_1d.cc +++ b/tests/grid/periodicity_1d.cc @@ -20,6 +20,7 @@ #include #include +#include #include #include "../tests.h" @@ -50,8 +51,10 @@ check() deallog << periodic_faces[i].cell[0]->index() << ' ' << periodic_faces[i].cell[1]->index() << ' ' << periodic_faces[i].face_idx[0] << ' ' - << periodic_faces[i].face_idx[1] << ' ' - << periodic_faces[i].orientation << std::endl; + << periodic_faces[i].face_idx[1] + << ' ' + // maintain old output by converting to a a bitset + << std::bitset<3>(periodic_faces[i].orientation) << std::endl; } } diff --git a/tests/mpi/periodicity_04.cc b/tests/mpi/periodicity_04.cc index e07bebfd03..33b84467e7 100644 --- a/tests/mpi/periodicity_04.cc +++ b/tests/mpi/periodicity_04.cc @@ -261,10 +261,10 @@ check(const unsigned int orientation, bool reverse) using CellFace = std::pair::cell_iterator, unsigned int>; - const typename std::map>> + const typename std::map> &face_map = triangulation.get_periodic_face_map(); typename std::map>>::const_iterator it; + std::pair>::const_iterator it; int sum_of_pairs_local = face_map.size(); int sum_of_pairs_global; MPI_Allreduce(&sum_of_pairs_local, @@ -295,7 +295,7 @@ check(const unsigned int orientation, bool reverse) { std::cout << "face_center_1: " << face_center_1 << std::endl; std::cout << "face_center_2: " << face_center_2 << std::endl; - typename std::map>>:: + typename std::map>:: const_iterator it; for (it = triangulation.get_periodic_face_map().begin(); it != triangulation.get_periodic_face_map().end(); -- 2.39.5