diff --git a/cpp/dolfinx/mesh/Topology.cpp b/cpp/dolfinx/mesh/Topology.cpp index 97fb34e96b..573fb122db 100644 --- a/cpp/dolfinx/mesh/Topology.cpp +++ b/cpp/dolfinx/mesh/Topology.cpp @@ -886,6 +886,20 @@ const std::vector& Topology::get_facet_permutations() const return _facet_permutations; } //----------------------------------------------------------------------------- +const std::vector& Topology::get_ridge_permutations() const +{ + if (auto i_map = this->index_map(this->dim() - 2); + !i_map + or (_ridge_permutations.empty() + and i_map->size_local() + i_map->num_ghosts() > 0)) + { + throw std::runtime_error( + "create_entity_permutations must be called before using this data."); + } + + return _ridge_permutations; +} +//----------------------------------------------------------------------------- const std::vector& Topology::interprocess_facets(int index) const { if (_interprocess_facets.empty()) @@ -1009,9 +1023,10 @@ void Topology::create_entity_permutations() for (int d = 0; d < tdim; ++d) create_entities(d); - auto [facet_permutations, cell_permutations] + auto [ridge_permutations, facet_permutations, cell_permutations] = compute_entity_permutations(*this); _facet_permutations = std::move(facet_permutations); + _ridge_permutations = std::move(ridge_permutations); _cell_permutations = std::move(cell_permutations); } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/mesh/Topology.h b/cpp/dolfinx/mesh/Topology.h index 681f88980e..12889fdf8d 100644 --- a/cpp/dolfinx/mesh/Topology.h +++ b/cpp/dolfinx/mesh/Topology.h @@ -159,6 +159,21 @@ class Topology /// computed const std::vector& get_facet_permutations() const; + /// @brief Get the numbers that encode the number of permutations to + /// apply to ridges. + /// + /// The permutations are encoded so that: + /// + /// @todo Matthew needs to document this. + /// + /// The data is stored in a flattened 2D array, so that `data[cell_index * + /// ridges_per_cell + ridge_index]` contains the ridge with index + /// `ridge_index` of the cell with index `cell_index`. + /// @return The encoded permutation info + /// @note An exception is raised if the permutations have not been + /// computed + const std::vector& get_ridge_permutations() const; + /// @brief Get the types of cells in the topology /// @return The cell types std::vector cell_types() const; @@ -236,6 +251,11 @@ class Topology // celln_0, celln_1, ,celln_2,] std::vector _facet_permutations; + // The ridge permutations (local ridge, cell) + // [cell0_0, cell0_1, ,cell0_2, cell1_0, cell1_1, ,cell1_2, ..., + // celln_0, celln_1, ,celln_2,] + std::vector _ridge_permutations; + // Cell permutation info. See the documentation for // get_cell_permutation_info for documentation of how this is encoded. std::vector _cell_permutations; diff --git a/cpp/dolfinx/mesh/permutationcomputation.cpp b/cpp/dolfinx/mesh/permutationcomputation.cpp index c64d4cc800..21a04080f8 100644 --- a/cpp/dolfinx/mesh/permutationcomputation.cpp +++ b/cpp/dolfinx/mesh/permutationcomputation.cpp @@ -284,17 +284,22 @@ compute_face_permutations(const mesh::Topology& topology) } // namespace //----------------------------------------------------------------------------- -std::pair, std::vector> +std::tuple, std::vector, + std::vector> mesh::compute_entity_permutations(const mesh::Topology& topology) { + common::Timer t_perm("Compute entity permutations"); const int tdim = topology.dim(); CellType cell_type = topology.cell_type(); const std::int32_t num_cells = topology.connectivity(tdim, 0)->num_nodes(); const int facets_per_cell = cell_num_entities(cell_type, tdim - 1); + const int ridges_per_cell = cell_num_entities(cell_type, tdim - 2); std::vector cell_permutation_info(num_cells, 0); std::vector facet_permutations(num_cells * facets_per_cell); + std::vector ridge_permutations(num_cells * ridges_per_cell); + std::int32_t used_bits = 0; if (tdim > 2) { @@ -335,9 +340,18 @@ mesh::compute_entity_permutations(const mesh::Topology& topology) facet_permutations[c * facets_per_cell + i] = edge_perm[c][i]; } } + if (tdim == 3) + { + for (int c = 0; c < num_cells; ++c) + { + for (int i = 0; i < ridges_per_cell; ++i) + ridge_permutations[c * ridges_per_cell + i] = edge_perm[c][i]; + } + } } assert(used_bits < _BITSETSIZE); - return {std::move(facet_permutations), std::move(cell_permutation_info)}; + return {std::move(ridge_permutations), std::move(facet_permutations), + std::move(cell_permutation_info)}; } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/mesh/permutationcomputation.h b/cpp/dolfinx/mesh/permutationcomputation.h index 1c38d09d2c..e1ccfedae4 100644 --- a/cpp/dolfinx/mesh/permutationcomputation.h +++ b/cpp/dolfinx/mesh/permutationcomputation.h @@ -66,8 +66,9 @@ class Topology; /// This data is used to correct the direction of vector function /// on permuted facets. /// -/// @return Facet permutation and cells permutations -std::pair, std::vector> +/// @param[in] topology The mesh topology. +std::tuple, std::vector, + std::vector> compute_entity_permutations(const Topology& topology); } // namespace dolfinx::mesh diff --git a/python/dolfinx/mesh.py b/python/dolfinx/mesh.py index a9ffe8d8ca..c9ec3606ec 100644 --- a/python/dolfinx/mesh.py +++ b/python/dolfinx/mesh.py @@ -173,6 +173,21 @@ def get_facet_permutations(self) -> npt.NDArray[np.uint8]: """ return self._cpp_object.get_facet_permutations() + def get_ridge_permutations(self) -> npt.NDArray[np.uint8]: + """Get the permutation integer to apply to ridges. + + The bits of each integer describes the number of reflections + that has to be applied to a ridge to map between a + ridge in the mesh (relative to a cell) and the corresponding + ridge on the reference element. The data has the shape + ``(num_cells, num_ridges_per_cell)``, flattened row-wise. The + number of cells include potential ghost cells. + + Note: + The data can be unpacked with ``numpy.unpack_bits``. + """ + return self._cpp_object.get_ridge_permutations() + def index_map(self, dim: int) -> _cpp.common.IndexMap: """Get the IndexMap that describes the parallel distribution of the mesh entities. diff --git a/python/dolfinx/wrappers/mesh.cpp b/python/dolfinx/wrappers/mesh.cpp index d74ded0c62..6483ef8dc7 100644 --- a/python/dolfinx/wrappers/mesh.cpp +++ b/python/dolfinx/wrappers/mesh.cpp @@ -668,6 +668,15 @@ void mesh(nb::module_& m) {p.size()}); }, nb::rv_policy::reference_internal) + .def( + "get_ridge_permutations", + [](const dolfinx::mesh::Topology& self) + { + const std::vector& p = self.get_ridge_permutations(); + return nb::ndarray(p.data(), + {p.size()}); + }, + nb::rv_policy::reference_internal) .def( "get_cell_permutation_info", [](const dolfinx::mesh::Topology& self)