diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ac68ebc406..e83999df91 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -164,6 +164,7 @@ target_sources( T8 PRIVATE t8_cmesh/t8_cmesh_internal/t8_cmesh_partition.cxx t8_cmesh/t8_cmesh_internal/t8_cmesh_stash.cxx t8_cmesh/t8_cmesh_internal/t8_cmesh_trees.cxx + t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_conn_tree_to_vertex.cxx t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_conn_vertex_to_tree.cxx t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_connectivity.cxx diff --git a/src/t8_cmesh/t8_cmesh.cxx b/src/t8_cmesh/t8_cmesh.cxx index efe0bd6ec8..d0ee721cd6 100644 --- a/src/t8_cmesh/t8_cmesh.cxx +++ b/src/t8_cmesh/t8_cmesh.cxx @@ -153,6 +153,21 @@ t8_cmesh_disable_negative_volume_check ([[maybe_unused]] t8_cmesh_t cmesh) #endif } +void +t8_cmesh_enable_tree_reordering (t8_cmesh_t cmesh, const int keep_reindex_map) +{ + T8_ASSERT (t8_cmesh_is_initialized (cmesh)); + cmesh->reindex_trees = 1; + cmesh->keep_reindex_map = keep_reindex_map; +} + +const t8_gloidx_t * +t8_cmesh_get_tree_reindex_map (const t8_cmesh_t cmesh) +{ + T8_ASSERT (t8_cmesh_is_committed (cmesh)); + return cmesh->tree_reindex_map; +} + #if T8_ENABLE_DEBUG int t8_cmesh_validate_geometry (const t8_cmesh_t cmesh, const int check_for_negative_volume) @@ -240,7 +255,9 @@ t8_cmesh_init (t8_cmesh_t *pcmesh) #if T8_ENABLE_DEBUG cmesh->negative_volume_check = 1; #endif /* T8_ENABLE_DEBUG */ - + cmesh->reindex_trees = 0; + cmesh->keep_reindex_map = 0; + cmesh->tree_reindex_map = nullptr; T8_ASSERT (t8_cmesh_is_initialized (cmesh)); } @@ -1248,6 +1265,11 @@ t8_cmesh_reset (t8_cmesh_t *pcmesh) cmesh->set_partition_scheme->unref (); } + /* Release the optional tree-ID mapping owned by the coarse mesh. */ + if (cmesh->tree_reindex_map != nullptr) { + T8_FREE (cmesh->tree_reindex_map); + } + delete cmesh->vertex_connectivity; delete cmesh->boundary_condition_handler; @@ -2261,8 +2283,8 @@ t8_cmesh_uniform_bounds_from_partition (const t8_cmesh_t cmesh, const t8_gloidx_ #endif } /* End sending of end message */ - } /* End loop over processes */ - } /* if (num_pure_local_trees > 0) */ + } /* End loop over processes */ + } /* if (num_pure_local_trees > 0) */ if (this_proc_is_empty) { /* If this proc is empty, we only expect a start message. diff --git a/src/t8_cmesh/t8_cmesh.h b/src/t8_cmesh/t8_cmesh.h index b5b6599904..a7f9bc9703 100644 --- a/src/t8_cmesh/t8_cmesh.h +++ b/src/t8_cmesh/t8_cmesh.h @@ -132,6 +132,22 @@ t8_cmesh_stash_is_empty (const t8_cmesh_t cmesh); void t8_cmesh_disable_negative_volume_check (t8_cmesh_t cmesh); +/** Enable locality-based coarse-tree reordering during \ref t8_cmesh_commit. + * \param [in,out] cmesh An initialized, uncommitted coarse mesh. + * \param [in] keep_reindex_map Nonzero to retain an old-to-new tree-ID mapping after commit. + * \note Enable on every process participating in commit. Stash records are gathered on rank zero, + * reordered there, and broadcast. Partition ranges and offsets refer to the reordered tree IDs. + */ +void +t8_cmesh_enable_tree_reordering (t8_cmesh_t cmesh, int keep_reindex_map); + +/** Return the stored old-to-new tree-ID mapping. + * \param [in] cmesh A committed coarse mesh. + * \return Read-only array indexed by the old global tree ID, or nullptr if none was retained. + * The array has t8_cmesh_get_num_trees(cmesh) entries and remains owned by \a cmesh. + */ +const t8_gloidx_t * +t8_cmesh_get_tree_reindex_map (const t8_cmesh_t cmesh); #if T8_ENABLE_DEBUG /** Check the geometry of the mesh for validity, this means checking if trees and their geometries * are compatible and if they have negative volume. @@ -140,7 +156,6 @@ t8_cmesh_disable_negative_volume_check (t8_cmesh_t cmesh); * \return True if the geometry of the cmesh is valid. */ int - t8_cmesh_validate_geometry (const t8_cmesh_t cmesh, const int check_for_negative_volume); #endif diff --git a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_commit.cxx b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_commit.cxx index 6ad3b414bc..91c292de72 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_commit.cxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_commit.cxx @@ -36,6 +36,7 @@ #include #include #include +#include #include /** @@ -545,6 +546,10 @@ t8_cmesh_commit_from_stash (t8_cmesh_t cmesh, sc_MPI_Comm comm) } #endif + if (cmesh->reindex_trees) { + t8_cmesh_reindex_trees (cmesh, comm); + } + if (cmesh->set_partition) { /* partitioned commit */ t8_cmesh_commit_partitioned_new (cmesh, comm); diff --git a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx new file mode 100644 index 0000000000..2a21f5dea5 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -0,0 +1,531 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + Copyright (C) 2026 the developers + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + General Public License for more details. + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ +/** \file t8_cmesh_tree_reindex.cxx + * Implements SFC-based reindexing of coarse mesh trees from their geometric vertex data. + */ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +/** \brief Stores the coarse-tree data assigned to one element of the auxiliary forest. + * + * Entries in \a tree_centers and \a tree_ids correspond by index. + */ +struct t8_cmesh_tree_reindex_element_data +{ + std::vector tree_centers; /**< Coordinates of the coarse-tree centers assigned to this element. */ + std::vector tree_ids; /**< Global coarse-tree IDs corresponding to \a tree_centers. */ +}; + +/** \brief Stores per-element data used while refining the auxiliary forest. + * + * The object owns an array with one entry for each local leaf element of the corresponding forest. + */ +struct t8_cmesh_tree_reindex_forest_data +{ + /** \brief Construct forest data for a given number of local leaf elements. + * \param [in] num_elements Number of local leaf elements for which element data is allocated. + * \param [in] tolerance Absolute tolerance used to determine whether a point lies inside a child. + */ + explicit t8_cmesh_tree_reindex_forest_data (const t8_locidx_t num_elements, const double tolerance) + : elements (T8_ALLOC (t8_cmesh_tree_reindex_element_data, num_elements)), num_elements (num_elements), + tolerance (tolerance), finished (true) + { + /* T8_ALLOC only reserves raw storage. Explicitly construct every C++ element-data object in that storage. */ + std::uninitialized_value_construct_n (elements, num_elements); + } + /** \brief Destroy all element data and release the owned array. */ + ~t8_cmesh_tree_reindex_forest_data () + { + /* Destroy the C++ objects before releasing the raw storage allocated by T8_ALLOC. */ + std::destroy_n (elements, num_elements); + T8_FREE (elements); + } + t8_cmesh_tree_reindex_element_data *elements; /**< Data for each forest-local leaf element. */ + t8_locidx_t num_elements; /**< Number of entries in \a elements. */ + double tolerance; /**< Geometric tolerance used while distributing points among refined children. */ + bool finished; /**< True if no local leaf element contains more than one coarse-tree center. */ +}; +/** + * \brief Refine an auxiliary forest element if it contains more than one coarse-tree center. + * + * This callback only requests refinement. It never requests coarsening or removal. + * + * \param [in] forest Forest to which the adapted elements will belong. Unused. + * \param [in] forest_from Source forest that is being adapted. + * \param [in] which_tree Local tree containing the considered element. + * \param [in] tree_class Element class of \a which_tree. + * \param [in] lelement_id Tree-local leaf index of the considered element in \a forest_from. + * \param [in] scheme Refinement scheme of the forest. + * \param [in] is_family Whether \a *elements* forms a family. Unused because no coarsening is performed. + * \param [in] num_elements Number of valid entries in \a *elements.* Unused. + * \param [in] elements Elements presented to the adaptation callback. Unused. + * \return 1 if the considered element contains more than one coarse-tree center, 0 otherwise. + */ +static int +t8_cmesh_tree_reindex_adapt ([[maybe_unused]] t8_forest_t forest, t8_forest_t forest_from, const t8_locidx_t which_tree, + const t8_eclass_t tree_class, const t8_locidx_t lelement_id, const t8_scheme *scheme, + [[maybe_unused]] const int is_family, [[maybe_unused]] const int num_elements, + [[maybe_unused]] t8_element_t *elements[]) +{ + const t8_cmesh_tree_reindex_forest_data *data + = static_cast (t8_forest_get_user_data (forest_from)); + T8_ASSERT (data != nullptr); + /* Convert the tree-local leaf index to the forest-local index used by the auxiliary data array. */ + const t8_locidx_t element_index = t8_forest_get_tree_element_offset (forest_from, which_tree) + lelement_id; + const t8_element_t *element = t8_forest_get_leaf_element_in_tree (forest_from, which_tree, lelement_id); + + /* Keep all centers in an element at the maximum refinement level for deterministic tie-breaking. */ + return data->elements[element_index].tree_centers.size () > 1 + && scheme->element_get_level (tree_class, element) < t8_forest_get_maxlevel (forest_from); +} + +/** + * \brief Transfer coarse-tree centers and IDs from an auxiliary forest to its adapted successor. + * + * Unchanged elements are copied directly. If an element was refined, each coarse-tree center is assigned to the + * first child element that contains it. The first-match rule prevents centers on child boundaries from being copied + * to more than one child. + * + * \param [in] forest_old Source forest before adaptation. + * \param [in] forest_new Destination forest after adaptation. + * \param [in] which_tree Local tree containing the replaced elements. + * \param [in] tree_class Element class of \a which_tree. + * \param [in] scheme Refinement scheme of the forest. + * \param [in] refine 0 if the element is unchanged, 1 if it was refined. + * \param [in] num_old_elements Number of replaced elements in \a *forest_old.* + * \param [in] first_old_element Tree-local index of the first replaced element in \a *forest_old.* + * \param [in] num_new_elements Number of replacement elements in \a *forest_new.* + * \param [in] first_new_element Tree-local index of the first replacement element in \a *forest_new.* + */ +static void +t8_cmesh_tree_reindex_replace (t8_forest_t forest_old, t8_forest_t forest_new, const t8_locidx_t which_tree, + const t8_eclass_t tree_class, const t8_scheme *scheme, const int refine, + [[maybe_unused]] const int num_old_elements, + [[maybe_unused]] const t8_locidx_t first_old_element, const int num_new_elements, + const t8_locidx_t first_new_element) +{ + const t8_cmesh_tree_reindex_forest_data *data_old + = static_cast (t8_forest_get_user_data (forest_old)); + t8_cmesh_tree_reindex_forest_data *data_new + = static_cast (t8_forest_get_user_data (forest_new)); + T8_ASSERT (data_old != nullptr); + T8_ASSERT (data_new != nullptr); + T8_ASSERT (refine == 0 || refine == 1); + const t8_locidx_t old_index = t8_forest_get_tree_element_offset (forest_old, which_tree) + first_old_element; + const t8_cmesh_tree_reindex_element_data &old_element_data = data_old->elements[old_index]; + T8_ASSERT (num_old_elements == 1); + T8_ASSERT (old_element_data.tree_centers.size () == old_element_data.tree_ids.size ()); + if (refine == 0) { + /* Unchanged elements have at most one center, except when the maximum refinement level was reached. */ + T8_ASSERT (num_new_elements == 1); + if (old_element_data.tree_centers.size () > 1) { + const t8_element_t *old_element = t8_forest_get_leaf_element_in_tree (forest_old, which_tree, first_old_element); + T8_ASSERT (scheme->element_get_level (tree_class, old_element) == t8_forest_get_maxlevel (forest_old)); + } + const t8_locidx_t new_index = t8_forest_get_tree_element_offset (forest_new, which_tree) + first_new_element; + data_new->elements[new_index] = old_element_data; + return; + } + /* Element was refined. So it should have more than one center. */ + T8_ASSERT (old_element_data.tree_centers.size () > 1); + /* t8_forest_element_points_inside expects the points as a flat xyz array. */ + std::vector point_coordinates; + point_coordinates.reserve (3 * old_element_data.tree_centers.size ()); + for (const t8_3D_vec &tree_center : old_element_data.tree_centers) { + point_coordinates.insert (point_coordinates.end (), tree_center.begin (), tree_center.end ()); + } + /* A center on a child boundary may be reported inside multiple children. Assign it only once. */ + std::vector point_was_copied (old_element_data.tree_centers.size (), 0); + for (int inew = 0; inew < num_new_elements; ++inew) { + const t8_locidx_t new_tree_leaf_index = first_new_element + inew; + const t8_locidx_t new_index = t8_forest_get_tree_element_offset (forest_new, which_tree) + new_tree_leaf_index; + const t8_element_t *new_element = t8_forest_get_leaf_element_in_tree (forest_new, which_tree, new_tree_leaf_index); + /* Determine which old tree centers lie in the current child element. */ + std::vector point_inside (old_element_data.tree_centers.size (), 0); + t8_forest_element_points_inside (forest_new, which_tree, new_element, point_coordinates.data (), + static_cast (old_element_data.tree_centers.size ()), point_inside.data (), + data_old->tolerance); + /* Copy each matching center together with the tree ID at the same vector index. */ + t8_cmesh_tree_reindex_element_data &new_element_data = data_new->elements[new_index]; + for (std::size_t ipoint = 0; ipoint < old_element_data.tree_centers.size (); ++ipoint) { + if (point_inside[ipoint] && !point_was_copied[ipoint]) { + new_element_data.tree_centers.push_back (old_element_data.tree_centers[ipoint]); + new_element_data.tree_ids.push_back (old_element_data.tree_ids[ipoint]); + point_was_copied[ipoint] = 1; + } + } + /* Request another pass only when an element contains multiple centers and can still be refined. */ + if (new_element_data.tree_centers.size () > 1 + && scheme->element_get_level (tree_class, new_element) < t8_forest_get_maxlevel (forest_new)) { + data_new->finished = false; + } + } + /* Every center of a refined parent must have been assigned to one of its children. */ + T8_ASSERT ( + std::all_of (point_was_copied.begin (), point_was_copied.end (), [] (const int copied) { return copied != 0; })); +} + +std::map +t8_cmesh_calculate_tree_reindexing_map (t8_cmesh_t cmesh, sc_MPI_Comm comm) +{ + /* Compute a geometric center for every tree and the bounding box of all vertices. The centers are inserted into + * a single-tree auxiliary forest, which is refined until each leaf contains at most one center or reaches its + * maximum level. Traversing the leaves in SFC order yields the global old-to-new tree-ID permutation. */ + T8_ASSERT (cmesh != nullptr); + T8_ASSERT (cmesh->stash != nullptr); + const t8_stash_t original_cmesh_stash = cmesh->stash; + std::map tree_to_center; + t8_3D_vec min_coordinates + = { std::numeric_limits::max (), std::numeric_limits::max (), std::numeric_limits::max () }; + t8_3D_vec max_coordinates = { std::numeric_limits::lowest (), std::numeric_limits::lowest (), + std::numeric_limits::lowest () }; + /* Compute each tree center and the bounding box of all stored tree vertices in one pass over the attributes. */ + for (size_t iattr = 0; iattr < original_cmesh_stash->attributes.elem_count; ++iattr) { + const t8_stash_attribute_struct_t *attribute + = static_cast (sc_array_index (&original_cmesh_stash->attributes, iattr)); + if (attribute->package_id != t8_get_package_id () || attribute->key != T8_CMESH_VERTICES_ATTRIBUTE_KEY) { + continue; + } + T8_ASSERT (attribute->attr_data != nullptr); + T8_ASSERT (attribute->attr_size % (3 * sizeof (double)) == 0); + const size_t num_vertices = attribute->attr_size / (3 * sizeof (double)); + T8_ASSERT (num_vertices > 0); + const t8_3D_vec *tree_vertices = static_cast (attribute->attr_data); + t8_3D_vec tree_center = {}; + for (size_t ivert = 0; ivert < num_vertices; ++ivert) { + const t8_3D_vec &vertex = tree_vertices[ivert]; + for (int idim = 0; idim < 3; ++idim) { + min_coordinates[idim] = std::min (min_coordinates[idim], vertex[idim]); + max_coordinates[idim] = std::max (max_coordinates[idim], vertex[idim]); + } + t8_axpy (vertex, tree_center, 1.0); + } + t8_ax (tree_center, 1.0 / static_cast (num_vertices)); + tree_to_center.emplace (attribute->id, tree_center); + } + SC_CHECK_ABORTF (tree_to_center.size () == original_cmesh_stash->classes.elem_count, + "SFC tree reindexing requires vertex attributes for every tree (%zu centers, %zu classes).", + tree_to_center.size (), original_cmesh_stash->classes.elem_count); + + t8_cmesh_t bbox_cmesh; + t8_cmesh_init (&bbox_cmesh); + /* Prevent committing the auxiliary cmesh from recursively invoking tree reindexing. */ + bbox_cmesh->reindex_trees = 0; + /* The axis-aligned geometry is described by the minimum and maximum corner only. */ + const std::array bbox_vertices = { min_coordinates[0], min_coordinates[1], min_coordinates[2], + max_coordinates[0], max_coordinates[1], max_coordinates[2] }; + const double dx = max_coordinates[0] - min_coordinates[0]; + const double dy = max_coordinates[1] - min_coordinates[1]; + const double dz = max_coordinates[2] - min_coordinates[2]; + const double max_extent = std::max ({ std::abs (dx), std::abs (dy), std::abs (dz) }); + const double tolerance = max_extent * T8_PRECISION_EPS; + /* Ignore numerically zero extents when selecting the element class of the bounding box. */ + const int active_dimensions = (std::abs (dx) > tolerance) + (std::abs (dy) > tolerance) + (std::abs (dz) > tolerance); + t8_eclass_t bbox_eclass; + switch (active_dimensions) { + case 3: + bbox_eclass = T8_ECLASS_HEX; + break; + case 2: + bbox_eclass = T8_ECLASS_QUAD; + break; + case 1: + bbox_eclass = T8_ECLASS_LINE; + break; + default: + SC_ABORT ("Bounding box has zero extent in all directions.\n"); + } + /* Build a single-tree auxiliary cmesh covering all original tree centers. */ + t8_cmesh_set_tree_class (bbox_cmesh, 0, bbox_eclass); + t8_cmesh_set_tree_vertices (bbox_cmesh, 0, bbox_vertices.data (), 2); + t8_cmesh_register_geometry (bbox_cmesh); + t8_cmesh_commit (bbox_cmesh, comm); + t8_forest_t bbox_forest = t8_forest_new_uniform (bbox_cmesh, t8_scheme_new_default (), 0, 0, comm); + auto *data = new t8_cmesh_tree_reindex_forest_data (t8_forest_get_local_num_leaf_elements (bbox_forest), tolerance); + T8_ASSERT (data->num_elements == 1); + /* Initially every coarse-tree center belongs to the single root element of the auxiliary forest. */ + data->elements[0].tree_centers.reserve (tree_to_center.size ()); + data->elements[0].tree_ids.reserve (tree_to_center.size ()); + for (const auto &[global_tree_id, tree_center] : tree_to_center) { + data->elements[0].tree_centers.push_back (tree_center); + data->elements[0].tree_ids.push_back (global_tree_id); + } + data->finished = data->elements[0].tree_centers.size () <= 1 || t8_forest_get_maxlevel (bbox_forest) == 0; + t8_forest_set_user_data (bbox_forest, data); + /* Refine repeatedly until every auxiliary leaf contains at most one coarse-tree center. */ + while (!data->finished) { + t8_forest_t adapted_forest; + t8_forest_init (&adapted_forest); + /* Refine exactly those leaves that still contain more than one center. */ + t8_forest_set_adapt (adapted_forest, bbox_forest, t8_cmesh_tree_reindex_adapt, 0); + t8_forest_ref (bbox_forest); + t8_forest_commit (adapted_forest); + /* Transfer the tree-center data from each replaced element to its successor elements. */ + auto *adapted_data + = new t8_cmesh_tree_reindex_forest_data (t8_forest_get_local_num_leaf_elements (adapted_forest), tolerance); + t8_forest_set_user_data (adapted_forest, adapted_data); + t8_forest_iterate_replace (adapted_forest, bbox_forest, t8_cmesh_tree_reindex_replace); + delete data; + t8_forest_unref (&bbox_forest); + bbox_forest = adapted_forest; + data = adapted_data; + } + std::map tree_reindex; + t8_gloidx_t new_tree_index = 0; + /* Forest leaves are stored in SFC order. Visiting them in this order defines the new tree indices. */ + const t8_locidx_t num_bbox_local_trees = t8_forest_get_num_local_trees (bbox_forest); + for (t8_locidx_t bbox_itree = 0; bbox_itree < num_bbox_local_trees; ++bbox_itree) { + const t8_locidx_t num_leaf_elements = t8_forest_get_tree_num_leaf_elements (bbox_forest, bbox_itree); + for (t8_locidx_t ielement = 0; ielement < num_leaf_elements; ++ielement) { + const t8_locidx_t element_index = t8_forest_get_tree_element_offset (bbox_forest, bbox_itree) + ielement; + const t8_cmesh_tree_reindex_element_data &element_data = data->elements[element_index]; + T8_ASSERT (element_data.tree_centers.size () == element_data.tree_ids.size ()); + if (element_data.tree_ids.empty ()) { + continue; + } + /* Several centers can remain together at the maximum refinement level (e.g. coincident centers). + * Preserve their original tree-ID order within this leaf instead of silently dropping them. */ + for (const t8_gloidx_t tree_id : element_data.tree_ids) { + tree_reindex.emplace (tree_id, new_tree_index); + ++new_tree_index; + } + } + } + SC_CHECK_ABORTF (tree_reindex.size () == tree_to_center.size (), + "SFC tree reindexing lost tree centers (%zu indexed, %zu expected).", tree_reindex.size (), + tree_to_center.size ()); + + delete data; + t8_forest_unref (&bbox_forest); + return tree_reindex; +} + +void +t8_cmesh_tree_perform_reindex_inplace (t8_stash_t &stash, const std::map &tree_reindex) +{ + T8_ASSERT (stash != nullptr); + /* Update the IDs attached to the stored tree classes. */ + for (size_t iclass = 0; iclass < stash->classes.elem_count; ++iclass) { + t8_stash_class_struct_t *tree_class + = static_cast (sc_array_index (&stash->classes, iclass)); + tree_class->id = tree_reindex.at (tree_class->id); + } + /* Attributes reference trees by ID as well and must be updated consistently. */ + for (size_t iattr = 0; iattr < stash->attributes.elem_count; ++iattr) { + t8_stash_attribute_struct_t *attribute + = static_cast (sc_array_index (&stash->attributes, iattr)); + attribute->id = tree_reindex.at (attribute->id); + } + /* Reindex both endpoints of every face connection. */ + for (size_t iface = 0; iface < stash->joinfaces.elem_count; ++iface) { + t8_stash_joinface_struct_t *join + = static_cast (sc_array_index (&stash->joinfaces, iface)); + join->id1 = tree_reindex.at (join->id1); + join->id2 = tree_reindex.at (join->id2); + /* Keep join faces in the canonical tree-ID order expected by the stash routines. */ + if (join->id1 > join->id2) { + std::swap (join->id1, join->id2); + std::swap (join->face1, join->face2); + } + } + /* Reindexing changes the key order of the stash arrays, so restore their canonical ordering. */ + t8_stash_class_sort (stash); + t8_stash_joinface_sort (stash); + t8_stash_attribute_sort (stash); +} + +/** Gather the stash records and attribute data on rank zero, removing replicated entries. + * \param [in,out] cmesh Coarse mesh whose stash is replaced by the complete stash on rank zero. + * \param [in] comm Communicator on which the stash records are gathered. + */ +static void +t8_cmesh_tree_reindex_gather_stash (t8_cmesh_t cmesh, sc_MPI_Comm comm) +{ + t8_stash_t gathered_stash; + t8_stash_init (&gathered_stash); + sc_array_t attribute_data; + sc_array_t gathered_data; + sc_array_init (&attribute_data, sizeof (char)); + sc_array_init (&gathered_data, sizeof (char)); + /* Pack attribute payloads separately because their pointers cannot be used on another process. */ + for (size_t iattr = 0; iattr < cmesh->stash->attributes.elem_count; ++iattr) { + const auto *attribute + = static_cast (sc_array_index (&cmesh->stash->attributes, iattr)); + if (attribute->attr_size > 0) { + memcpy (sc_array_push_count (&attribute_data, attribute->attr_size), attribute->attr_data, attribute->attr_size); + } + } + + sc_array_t *local_arrays[] + = { &cmesh->stash->classes, &cmesh->stash->joinfaces, &cmesh->stash->attributes, &attribute_data }; + sc_array_t *gathered_arrays[] + = { &gathered_stash->classes, &gathered_stash->joinfaces, &gathered_stash->attributes, &gathered_data }; + std::vector recvcounts (cmesh->mpisize); + std::vector displacements (cmesh->mpisize); + /* Gather byte counts first to size the receive buffers for inputs of different lengths. */ + for (int iarray = 0; iarray < 4; ++iarray) { + const sc_array_t *local_array = local_arrays[iarray]; + sc_array_t *gathered_array = gathered_arrays[iarray]; + const size_t num_bytes = local_array->elem_count * local_array->elem_size; + SC_CHECK_ABORT (num_bytes <= static_cast (std::numeric_limits::max ()), + "Tree reindexing gather exceeds the MPI count limit."); + int sendcount = static_cast (num_bytes); + SC_CHECK_MPI (sc_MPI_Gather (&sendcount, 1, sc_MPI_INT, recvcounts.data (), 1, sc_MPI_INT, 0, comm)); + if (cmesh->mpirank == 0) { + int total_bytes = 0; + for (int iproc = 0; iproc < cmesh->mpisize; ++iproc) { + displacements[iproc] = total_bytes; + SC_CHECK_ABORT (recvcounts[iproc] <= std::numeric_limits::max () - total_bytes, + "Tree reindexing gather exceeds the MPI count limit."); + total_bytes += recvcounts[iproc]; + } + sc_array_resize (gathered_array, total_bytes / gathered_array->elem_size); + } + SC_CHECK_MPI (sc_MPI_Gatherv (local_array->array, sendcount, sc_MPI_BYTE, gathered_array->array, recvcounts.data (), + displacements.data (), sc_MPI_BYTE, 0, comm)); + } + + if (cmesh->mpirank == 0) { + /* Replicated stashes and ghost records can contain duplicates; keep one copy of each entry. */ + t8_stash_destroy (&cmesh->stash); + t8_stash_init (&cmesh->stash); + std::map tree_classes; + for (size_t iclass = 0; iclass < gathered_stash->classes.elem_count; ++iclass) { + const auto *tree_class + = static_cast (sc_array_index (&gathered_stash->classes, iclass)); + const auto [entry, inserted] = tree_classes.emplace (tree_class->id, tree_class->eclass); + SC_CHECK_ABORT (entry->second == tree_class->eclass, "Conflicting tree classes in reindexing input."); + if (inserted) { + t8_stash_add_class (cmesh->stash, tree_class->id, tree_class->eclass); + } + } + std::set> joins; + for (size_t ijoin = 0; ijoin < gathered_stash->joinfaces.elem_count; ++ijoin) { + const auto *join + = static_cast (sc_array_index (&gathered_stash->joinfaces, ijoin)); + if (joins.emplace (join->id1, join->id2, join->face1, join->face2, join->orientation).second) { + t8_stash_add_facejoin (cmesh->stash, join->id1, join->id2, join->face1, join->face2, join->orientation); + } + } + std::map, const t8_stash_attribute_struct_t *> attributes; + /* Attribute metadata and payloads were gathered in the same process order. */ + size_t data_offset = 0; + for (size_t iattr = 0; iattr < gathered_stash->attributes.elem_count; ++iattr) { + auto *attribute + = static_cast (sc_array_index (&gathered_stash->attributes, iattr)); + /* The received pointer belongs to the sending process. Copy the separately gathered payload instead. */ + attribute->is_owned = 0; + attribute->attr_data = gathered_data.array + data_offset; + const auto [entry, inserted] + = attributes.emplace (std::make_tuple (attribute->id, attribute->package_id, attribute->key), attribute); + if (inserted) { + t8_stash_add_attribute (cmesh->stash, attribute->id, attribute->package_id, attribute->key, + attribute->attr_size, attribute->attr_data, 1); + } + else { + const t8_stash_attribute_struct_t *previous = entry->second; + SC_CHECK_ABORT (previous->attr_size == attribute->attr_size + && (attribute->attr_size == 0 + || memcmp (previous->attr_data, attribute->attr_data, attribute->attr_size) == 0), + "Conflicting tree attributes in reindexing input."); + } + data_offset += attribute->attr_size; + } + T8_ASSERT (data_offset == gathered_data.elem_count); + } + t8_stash_destroy (&gathered_stash); + sc_array_reset (&attribute_data); + sc_array_reset (&gathered_data); +} + +void +t8_cmesh_reindex_trees (t8_cmesh_t cmesh, sc_MPI_Comm comm) +{ + T8_ASSERT (cmesh != nullptr); + T8_ASSERT (cmesh->stash != nullptr); + T8_ASSERT (cmesh->tree_reindex_map == nullptr); + + /* All processes must join the map broadcast if any process requests retention. */ + int local_keep_reindex_map = cmesh->keep_reindex_map != 0; + int keep_reindex_map = 0; + SC_CHECK_MPI (sc_MPI_Allreduce (&local_keep_reindex_map, &keep_reindex_map, 1, sc_MPI_INT, sc_MPI_MAX, comm)); + /* Empty input stashes also participate so rank zero receives the complete mesh. */ + if (cmesh->mpisize > 1) { + t8_cmesh_tree_reindex_gather_stash (cmesh, comm); + } + size_t elem_counts[3] = {}; + /* Calculate and apply the global permutation only on rank zero. */ + if (cmesh->mpirank == 0) { + const size_t num_trees = cmesh->stash->classes.elem_count; + if (num_trees > 0) { + const std::map reindexing_map = t8_cmesh_calculate_tree_reindexing_map (cmesh); + if (keep_reindex_map) { + SC_CHECK_ABORT (reindexing_map.size () == num_trees, "Incomplete SFC reindexing map."); + cmesh->tree_reindex_map = T8_ALLOC (t8_gloidx_t, num_trees); + /* Store a direct lookup from the original global tree ID to its new ID for user data. */ + for (const auto &[old_tree_id, new_tree_id] : reindexing_map) { + SC_CHECK_ABORTF (old_tree_id >= 0 && static_cast (old_tree_id) < num_trees, + "Invalid old tree ID %" T8_GLOIDX_FORMAT " in SFC reindexing map.", old_tree_id); + cmesh->tree_reindex_map[old_tree_id] = new_tree_id; + } + } + t8_cmesh_tree_perform_reindex_inplace (cmesh->stash, reindexing_map); + } + /* The stash broadcast expects counts in attribute, class, and join order. */ + elem_counts[0] = cmesh->stash->attributes.elem_count; + elem_counts[1] = num_trees; + elem_counts[2] = cmesh->stash->joinfaces.elem_count; + } + if (cmesh->mpisize == 1) { + return; + } + SC_CHECK_MPI (sc_MPI_Bcast (elem_counts, sizeof (elem_counts), sc_MPI_BYTE, 0, comm)); + /* Release the original attributes before receiving the reordered stash. */ + if (cmesh->mpirank != 0) { + t8_stash_destroy (&cmesh->stash); + t8_stash_init (&cmesh->stash); + } + t8_stash_bcast (cmesh->stash, 0, comm, elem_counts); + if (keep_reindex_map && elem_counts[1] > 0) { + SC_CHECK_ABORT (elem_counts[1] <= static_cast (std::numeric_limits::max ()), + "Tree reindexing map exceeds the MPI count limit."); + if (cmesh->mpirank != 0) { + cmesh->tree_reindex_map = T8_ALLOC (t8_gloidx_t, elem_counts[1]); + } + SC_CHECK_MPI (sc_MPI_Bcast (cmesh->tree_reindex_map, static_cast (elem_counts[1]), T8_MPI_GLOIDX, 0, comm)); + /* Discard the temporary map on processes that did not request retention. */ + if (!cmesh->keep_reindex_map) { + T8_FREE (cmesh->tree_reindex_map); + cmesh->tree_reindex_map = nullptr; + } + } +} diff --git a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.hxx b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.hxx new file mode 100644 index 0000000000..5e978aa59b --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.hxx @@ -0,0 +1,76 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2026 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** \file t8_cmesh_tree_reindex.hxx + * Internal declarations for reindexing coarse trees by geometric locality. + */ + +#ifndef T8_CMESH_TREE_REINDEX_HXX +#define T8_CMESH_TREE_REINDEX_HXX + +#include +#include + +#include + +/** Calculate the old-to-new coarse-tree ID mapping from the auxiliary forest's SFC order. + * + * The tree centers are computed from the vertex attributes of all coarse trees. + * The auxiliary forest is adapted until its leaves isolate the centers or reach + * the maximum refinement level. Multiple centers in one terminal leaf are ordered + * by their original global tree IDs. + * + * \param [in] cmesh Uncommitted coarse mesh with a complete replicated stash. + * \param [in] comm Communicator used for the auxiliary forest. Defaults to MPI self. + * \return Map from original global tree ID to its new SFC-based ID. + * + * \note This function only calculates the permutation and does not modify \a cmesh. + */ +std::map +t8_cmesh_calculate_tree_reindexing_map (t8_cmesh_t cmesh, sc_MPI_Comm comm = sc_MPI_COMM_SELF); + +/** Apply an already calculated old-to-new tree-ID permutation to a stash. + * Changes tree classes, attributes, and both endpoints of every face join, then + * restores the canonical ordering of the stash entries. + * + * \param [in,out] stash Stash containing the tree records to reindex. + * \param [in] tree_reindex Map from old global IDs to new global IDs. + */ +void +t8_cmesh_tree_perform_reindex_inplace (t8_stash_t &stash, const std::map &tree_reindex); + +/** Calculate and apply SFC-based coarse-tree reindexing. + * + * Gathers the stash on rank zero, calls \ref t8_cmesh_calculate_tree_reindexing_map and + * \ref t8_cmesh_tree_perform_reindex_inplace there, then broadcasts the reordered stash. + * If requested with + * \ref t8_cmesh_enable_tree_reordering, the original-to-new mapping is + * retained in the coarse mesh for retrieval after commit. + * + * \param [in,out] cmesh Uncommitted coarse mesh with replicated or distributed stash records. + * \param [in] comm Communicator on which all processes participate in reindexing. + * \note For partitioned coarse meshes, partition ranges and offsets refer to the reordered tree IDs. + */ +void +t8_cmesh_reindex_trees (t8_cmesh_t cmesh, sc_MPI_Comm comm); + +#endif /* T8_CMESH_TREE_REINDEX_HXX */ diff --git a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_types.h b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_types.h index 5f143865ff..292888058b 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_types.h +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_types.h @@ -101,6 +101,10 @@ typedef struct t8_cmesh int face_knowledge; /**< If partitioned the level of face knowledge that is expected. */ + int reindex_trees; /**< If nonzero the trees will be indexed according to their geometric locality */ + int keep_reindex_map; /**< Retain the old-to-new tree-ID mapping after commit if nonzero. */ + t8_gloidx_t *tree_reindex_map; /**< Retained mapping indexed by original global tree ID, or nullptr. */ + const t8_scheme_c *set_partition_scheme; /**< If the cmesh is to be partitioned according to a uniform level, the scheme that describes the refinement pattern. */ int8_t set_partition_level; /**< Non-negative if the cmesh should be partitioned from an already existing cmesh diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3365b6ccdd..a7e1995fa2 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -136,6 +136,8 @@ add_t8_cpp_test( NAME t8_gtest_cmesh_add_attributes_when_derive_parallel SOUR add_t8_cpp_test( NAME t8_gtest_cmesh_vertex_conn_tree_to_vertex_parallel SOURCES t8_cmesh/t8_gtest_cmesh_vertex_conn_tree_to_vertex.cxx ) add_t8_cpp_test( NAME t8_gtest_cmesh_vertex_conn_vertex_to_tree_parallel SOURCES t8_cmesh/t8_gtest_cmesh_vertex_conn_vertex_to_tree.cxx ) add_t8_cpp_test( NAME t8_gtest_cmesh_vertex_conn_serial SOURCES t8_cmesh/t8_gtest_cmesh_vertex_conn.cxx ) +add_t8_cpp_test( NAME t8_gtest_cmesh_tree_reindex_serial SOURCES t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx) +add_t8_cpp_test( NAME t8_gtest_cmesh_tree_reindex_parallel SOURCES t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx ) add_t8_cpp_test( NAME t8_gtest_multiple_attributes_parallel SOURCES t8_cmesh/t8_gtest_multiple_attributes.cxx ) add_t8_cpp_test( NAME t8_gtest_attribute_gloidx_array_serial SOURCES t8_cmesh/t8_gtest_attribute_gloidx_array.cxx ) add_t8_cpp_test( NAME t8_gtest_cmesh_bounding_box_serial SOURCES t8_cmesh/t8_gtest_cmesh_bounding_box.cxx ) diff --git a/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx b/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx new file mode 100644 index 0000000000..a04ed4c201 --- /dev/null +++ b/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx @@ -0,0 +1,495 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2026 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** \file t8_gtest_cmesh_tree_reindex.cxx + * Test opt-in tree reordering and the optional old-to-new tree-ID map. + */ + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +/* -1 keeps the default (disabled), 0 enables reordering without retaining the map, + * and 1 enables reordering and retains the map. */ +struct t8_test_cmesh_tree_reindex: public testing::TestWithParam +{ + protected: + void + SetUp () override + { + t8_cmesh_init (&cmesh); + t8_cmesh_init (&control_cmesh); + if (GetParam () >= 0) { + t8_cmesh_enable_tree_reordering (cmesh, GetParam ()); + } + } + + void + TearDown () override + { + t8_cmesh_unref (&cmesh); + t8_cmesh_unref (&control_cmesh); + } + + static void + build_lines (t8_cmesh_t mesh, const std::vector> &endpoints) + { + for (t8_gloidx_t tree_id = 0; tree_id < static_cast (endpoints.size ()); ++tree_id) { + const double vertices[6] = { endpoints[tree_id][0], 0, 0, endpoints[tree_id][1], 0, 0 }; + t8_cmesh_set_tree_class (mesh, tree_id, T8_ECLASS_LINE); + t8_cmesh_set_tree_vertices (mesh, tree_id, vertices, 2); + t8_cmesh_set_attribute (mesh, tree_id, t8_testsuite_get_package_id (), 0, &tree_id, sizeof (tree_id), 0); + + /* Join touching line segments, even when their original tree IDs are out of order. */ + for (t8_gloidx_t neighbor = 0; neighbor < static_cast (endpoints.size ()); ++neighbor) { + if (endpoints[tree_id][1] == endpoints[neighbor][0]) { + t8_cmesh_set_join (mesh, tree_id, neighbor, 1, 0, 0); + } + } + } + t8_cmesh_register_geometry (mesh); + t8_cmesh_commit (mesh, sc_MPI_COMM_SELF); + } + + void + check_mesh (const std::vector &expected_reindex) + { + ASSERT_TRUE (t8_cmesh_is_committed (cmesh)); + ASSERT_TRUE (t8_cmesh_is_committed (control_cmesh)); + const t8_gloidx_t num_trees = t8_cmesh_get_num_trees (control_cmesh); + ASSERT_EQ (expected_reindex.size (), static_cast (num_trees)); + ASSERT_EQ (t8_cmesh_get_num_trees (cmesh), num_trees); + ASSERT_EQ (t8_cmesh_get_num_local_trees (cmesh), num_trees); + EXPECT_EQ (t8_cmesh_get_tree_reindex_map (control_cmesh), nullptr); + + const t8_gloidx_t *reindex_map = t8_cmesh_get_tree_reindex_map (cmesh); + if (GetParam () == 1) { + ASSERT_NE (reindex_map, nullptr); + } + else { + EXPECT_EQ (reindex_map, nullptr); + } + + /* Check the permutation before using it to look up trees or neighbors. */ + std::vector seen (num_trees, false); + for (const t8_gloidx_t new_tree_id : expected_reindex) { + ASSERT_GE (new_tree_id, 0); + ASSERT_LT (new_tree_id, num_trees); + ASSERT_FALSE (seen[new_tree_id]); + seen[new_tree_id] = true; + } + + for (t8_gloidx_t old_tree_id = 0; old_tree_id < num_trees; ++old_tree_id) { + const t8_gloidx_t new_tree_id = expected_reindex[old_tree_id]; + SCOPED_TRACE (testing::Message () << "old tree " << old_tree_id << ", new tree " << new_tree_id); + if (reindex_map != nullptr) { + EXPECT_EQ (reindex_map[old_tree_id], new_tree_id); + } + const t8_locidx_t old_local_id = t8_cmesh_get_local_id (control_cmesh, old_tree_id); + const t8_locidx_t new_local_id = t8_cmesh_get_local_id (cmesh, new_tree_id); + ASSERT_GE (old_local_id, 0); + ASSERT_GE (new_local_id, 0); + EXPECT_EQ (t8_cmesh_get_global_id (cmesh, new_local_id), new_tree_id); + const t8_eclass_t eclass = t8_cmesh_get_tree_class (control_cmesh, old_local_id); + ASSERT_EQ (t8_cmesh_get_tree_class (cmesh, new_local_id), eclass); + + const double *old_vertices = t8_cmesh_get_tree_vertices (control_cmesh, old_local_id); + const double *new_vertices = t8_cmesh_get_tree_vertices (cmesh, new_local_id); + ASSERT_NE (old_vertices, nullptr); + ASSERT_NE (new_vertices, nullptr); + for (int icoord = 0; icoord < 3 * t8_eclass_num_vertices[eclass]; ++icoord) { + EXPECT_EQ (new_vertices[icoord], old_vertices[icoord]); + } + + /* Check every directed face, including boundaries and reverse joins. */ + for (int face = 0; face < t8_eclass_num_faces[eclass]; ++face) { + SCOPED_TRACE (testing::Message () << "face " << face); + int old_dual_face = -1, old_orientation = -1; + int new_dual_face = -1, new_orientation = -1; + const t8_locidx_t old_neighbor + = t8_cmesh_get_face_neighbor (control_cmesh, old_local_id, face, &old_dual_face, &old_orientation); + const t8_locidx_t new_neighbor + = t8_cmesh_get_face_neighbor (cmesh, new_local_id, face, &new_dual_face, &new_orientation); + if (old_neighbor < 0) { + EXPECT_LT (new_neighbor, 0); + } + else { + ASSERT_GE (new_neighbor, 0); + const t8_gloidx_t old_neighbor_global = t8_cmesh_get_global_id (control_cmesh, old_neighbor); + EXPECT_EQ (t8_cmesh_get_global_id (cmesh, new_neighbor), expected_reindex[old_neighbor_global]); + EXPECT_EQ (new_dual_face, old_dual_face); + EXPECT_EQ (new_orientation, old_orientation); + } + } + } + + /* Retrieving the map must remain possible after querying the committed mesh. */ + const t8_gloidx_t *retrieved_map = t8_cmesh_get_tree_reindex_map (cmesh); + ASSERT_EQ (retrieved_map, reindex_map); + if (retrieved_map != nullptr) { + for (t8_gloidx_t old_tree_id = 0; old_tree_id < num_trees; ++old_tree_id) { + EXPECT_EQ (retrieved_map[old_tree_id], expected_reindex[old_tree_id]); + } + } + } + + void + check_lines (const std::vector> &endpoints, std::vector expected_reindex) + { + build_lines (control_cmesh, endpoints); + build_lines (cmesh, endpoints); + if (GetParam () < 0) { + std::iota (expected_reindex.begin (), expected_reindex.end (), 0); + } + ASSERT_NO_FATAL_FAILURE (check_mesh (expected_reindex)); + + /* User attributes must stay attached to the original tree. */ + for (t8_gloidx_t old_tree_id = 0; old_tree_id < static_cast (endpoints.size ()); ++old_tree_id) { + const t8_locidx_t new_local_id = t8_cmesh_get_local_id (cmesh, expected_reindex[old_tree_id]); + const auto *original_id = static_cast ( + t8_cmesh_get_attribute (cmesh, t8_testsuite_get_package_id (), 0, new_local_id)); + ASSERT_NE (original_id, nullptr); + EXPECT_EQ (*original_id, old_tree_id); + } + } + + /** Compare a partitioned hypercube against the same ordering on a replicated mesh. + * \param [in] eclass Element class used for both hypercubes. + * \param [in] periodic Whether to connect opposite domain faces. + */ + void + check_partitioned_hypercube (const t8_eclass_t eclass, const int periodic) + { + int mpirank; + int mpisize; + SC_CHECK_MPI (sc_MPI_Comm_rank (sc_MPI_COMM_WORLD, &mpirank)); + SC_CHECK_MPI (sc_MPI_Comm_size (sc_MPI_COMM_WORLD, &mpisize)); + + /* Exercise overlapping input records by supplying the complete stash on every rank. */ + if (GetParam () >= 0) { + t8_cmesh_enable_tree_reordering (control_cmesh, 1); + } + t8_cmesh_new_hypercube (&control_cmesh, eclass, sc_MPI_COMM_SELF, 0, 0, periodic); + t8_cmesh_new_hypercube (&cmesh, eclass, sc_MPI_COMM_WORLD, 0, 1, periodic); + ASSERT_TRUE (t8_cmesh_is_partitioned (cmesh)); + const t8_gloidx_t num_trees = t8_cmesh_get_num_trees (control_cmesh); + const t8_gloidx_t first_tree = mpirank * num_trees / mpisize; + const t8_locidx_t num_local = (mpirank + 1) * num_trees / mpisize - first_tree; + check_partitioned_mesh (first_tree, num_local); + } + + /** Compare local trees and ghosts with a replicated reference mesh. + * \param [in] first_tree Expected first local global tree ID. + * \param [in] num_local Expected number of local trees. + */ + void + check_partitioned_mesh (const t8_gloidx_t first_tree, const t8_locidx_t num_local) + { + int mpirank; + SC_CHECK_MPI (sc_MPI_Comm_rank (sc_MPI_COMM_WORLD, &mpirank)); + ASSERT_TRUE (t8_cmesh_is_partitioned (cmesh)); + const t8_gloidx_t num_trees = t8_cmesh_get_num_trees (control_cmesh); + ASSERT_EQ (t8_cmesh_get_num_trees (cmesh), num_trees); + ASSERT_EQ (t8_cmesh_get_first_treeid (cmesh), first_tree); + ASSERT_EQ (t8_cmesh_get_num_local_trees (cmesh), num_local); + + /* Check the full map on every rank, including ranks without local trees. */ + const t8_gloidx_t *expected_map = t8_cmesh_get_tree_reindex_map (control_cmesh); + const t8_gloidx_t *reindex_map = t8_cmesh_get_tree_reindex_map (cmesh); + if (GetParam () == 1) { + ASSERT_NE (expected_map, nullptr); + ASSERT_NE (reindex_map, nullptr); + } + else { + EXPECT_EQ (reindex_map, nullptr); + } + + std::vector expected_ghosts (num_trees, false); + for (t8_locidx_t itree = 0; itree < num_local + t8_cmesh_get_num_ghosts (cmesh); ++itree) { + const t8_gloidx_t global_tree_id = t8_cmesh_get_global_id (cmesh, itree); + ASSERT_GE (global_tree_id, 0); + ASSERT_LT (global_tree_id, num_trees); + SCOPED_TRACE (testing::Message () << "rank " << mpirank << ", tree " << global_tree_id); + const t8_eclass_t eclass = t8_cmesh_get_tree_class (control_cmesh, global_tree_id); + if (itree >= num_local) { + EXPECT_EQ (t8_cmesh_get_ghost_class (cmesh, itree - num_local), eclass); + EXPECT_TRUE (expected_ghosts[global_tree_id]); + expected_ghosts[global_tree_id] = false; + } + else { + EXPECT_EQ (t8_cmesh_get_tree_class (cmesh, itree), eclass); + } + const double *vertices = t8_cmesh_get_tree_vertices (cmesh, itree); + const double *expected_vertices = t8_cmesh_get_tree_vertices (control_cmesh, global_tree_id); + ASSERT_NE (vertices, nullptr); + ASSERT_NE (expected_vertices, nullptr); + for (int icoordinate = 0; icoordinate < 3 * t8_eclass_num_vertices[eclass]; ++icoordinate) { + EXPECT_EQ (vertices[icoordinate], expected_vertices[icoordinate]); + } + for (int iface = 0; iface < t8_eclass_num_faces[eclass]; ++iface) { + int dual_face = -1; + int orientation = -1; + int expected_dual_face = -1; + int expected_orientation = -1; + const t8_locidx_t reference_neighbor = t8_cmesh_get_face_neighbor (control_cmesh, global_tree_id, iface, + &expected_dual_face, &expected_orientation); + const t8_locidx_t neighbor = t8_cmesh_get_face_neighbor (cmesh, itree, iface, &dual_face, &orientation); + EXPECT_EQ (t8_cmesh_tree_face_is_boundary (cmesh, itree, iface), + t8_cmesh_tree_face_is_boundary (control_cmesh, global_tree_id, iface)); + if (itree < num_local && reference_neighbor >= 0 + && (reference_neighbor < first_tree || reference_neighbor >= first_tree + num_local)) { + expected_ghosts[reference_neighbor] = true; + } + const t8_locidx_t expected_neighbor + = reference_neighbor < 0 ? -1 : t8_cmesh_get_local_id (cmesh, reference_neighbor); + if (expected_neighbor < 0) { + EXPECT_LT (neighbor, 0); + } + else { + EXPECT_EQ (neighbor, expected_neighbor); + EXPECT_EQ (dual_face, expected_dual_face); + EXPECT_EQ (orientation, expected_orientation); + } + } + } + EXPECT_EQ (std::count (expected_ghosts.begin (), expected_ghosts.end (), true), 0); + ASSERT_EQ (t8_cmesh_get_tree_reindex_map (cmesh), reindex_map); + if (reindex_map != nullptr) { + ASSERT_NE (expected_map, nullptr); + for (t8_gloidx_t itree = 0; itree < num_trees; ++itree) { + EXPECT_EQ (reindex_map[itree], expected_map[itree]); + } + } + } + + /** Check distributed line records, including copied and borrowed variable-size attributes. + * \param [in] concentrate Whether only the last rank supplies and owns trees. + */ + void + check_distributed_lines (const bool concentrate) + { + int mpirank; + int mpisize; + SC_CHECK_MPI (sc_MPI_Comm_rank (sc_MPI_COMM_WORLD, &mpirank)); + SC_CHECK_MPI (sc_MPI_Comm_size (sc_MPI_COMM_WORLD, &mpisize)); + const std::vector> endpoints = { { 3, 4 }, { 0, 1 }, { 5, 6 }, { 1, 2 }, { 4, 5 }, { 2, 3 } }; + std::vector expected_reindex = { 3, 0, 5, 1, 4, 2 }; + const t8_gloidx_t num_trees = endpoints.size (); + const t8_gloidx_t first_tree = concentrate ? 0 : mpirank * num_trees / mpisize; + const t8_locidx_t num_local + = concentrate ? (mpirank == mpisize - 1 ? num_trees : 0) : (mpirank + 1) * num_trees / mpisize - first_tree; + if (GetParam () >= 0) { + t8_cmesh_enable_tree_reordering (control_cmesh, 1); + } + else { + std::iota (expected_reindex.begin (), expected_reindex.end (), 0); + } + build_lines (control_cmesh, endpoints); + + /* Keep only the original local trees and their face ghosts in each input stash. */ + std::vector supplied_trees (num_trees, false); + for (t8_gloidx_t itree = first_tree; itree < first_tree + num_local; ++itree) { + supplied_trees[itree] = true; + for (t8_gloidx_t ineighbor = 0; ineighbor < num_trees; ++ineighbor) { + if (endpoints[itree][0] == endpoints[ineighbor][1] || endpoints[itree][1] == endpoints[ineighbor][0]) { + supplied_trees[ineighbor] = true; + } + } + } + std::vector> attributes (num_trees); + t8_cmesh_set_dimension (cmesh, 1); + for (t8_gloidx_t itree = 0; itree < num_trees; ++itree) { + attributes[itree].resize (3 * itree + 1); + for (size_t ibyte = 0; ibyte < attributes[itree].size (); ++ibyte) { + attributes[itree][ibyte] = static_cast (17 * itree + ibyte); + } + if (supplied_trees[itree]) { + const double vertices[6] = { endpoints[itree][0], 0, 0, endpoints[itree][1], 0, 0 }; + t8_cmesh_set_tree_class (cmesh, itree, T8_ECLASS_LINE); + t8_cmesh_set_tree_vertices (cmesh, itree, vertices, 2); + t8_cmesh_set_attribute (cmesh, itree, t8_testsuite_get_package_id (), 0, &itree, sizeof (itree), 0); + t8_cmesh_set_attribute (cmesh, itree, t8_testsuite_get_package_id (), 1, attributes[itree].data (), + attributes[itree].size (), itree % 2); + } + /* Include outward ghost joins whose remote endpoint may be absent from the local stash. */ + for (t8_gloidx_t ineighbor = 0; ineighbor < num_trees; ++ineighbor) { + if (endpoints[itree][1] == endpoints[ineighbor][0] && (supplied_trees[itree] || supplied_trees[ineighbor])) { + t8_cmesh_set_join (cmesh, itree, ineighbor, 1, 0, 0); + } + } + } + t8_cmesh_set_partition_range (cmesh, 3, first_tree, first_tree + num_local - 1); + t8_cmesh_register_geometry (cmesh); + t8_cmesh_commit (cmesh, sc_MPI_COMM_WORLD); + ASSERT_NO_FATAL_FAILURE (check_partitioned_mesh (first_tree, num_local)); + + const t8_gloidx_t *reindex_map = t8_cmesh_get_tree_reindex_map (cmesh); + if (reindex_map != nullptr) { + for (t8_gloidx_t itree = 0; itree < num_trees; ++itree) { + EXPECT_EQ (reindex_map[itree], expected_reindex[itree]); + } + } + for (t8_locidx_t itree = 0; itree < num_local + t8_cmesh_get_num_ghosts (cmesh); ++itree) { + const t8_gloidx_t *original_id + = static_cast (t8_cmesh_get_attribute (cmesh, t8_testsuite_get_package_id (), 0, itree)); + ASSERT_NE (original_id, nullptr); + ASSERT_GE (*original_id, 0); + ASSERT_LT (*original_id, num_trees); + EXPECT_EQ (t8_cmesh_get_global_id (cmesh, itree), expected_reindex[*original_id]); + const unsigned char *attribute + = static_cast (t8_cmesh_get_attribute (cmesh, t8_testsuite_get_package_id (), 1, itree)); + ASSERT_NE (attribute, nullptr); + for (size_t ibyte = 0; ibyte < attributes[*original_id].size (); ++ibyte) { + EXPECT_EQ (attribute[ibyte], attributes[*original_id][ibyte]); + } + } + EXPECT_EQ (t8_cmesh_get_tree_reindex_map (cmesh), reindex_map); + } + + t8_cmesh_t cmesh = nullptr; + t8_cmesh_t control_cmesh = nullptr; +}; + +TEST_P (t8_test_cmesh_tree_reindex, scrambled_lines) +{ + /* The 1D SFC orders these segments from left to right. */ + check_lines ({ { 3, 4 }, { 0, 1 }, { 5, 6 }, { 1, 2 }, { 4, 5 }, { 2, 3 } }, { 3, 0, 5, 1, 4, 2 }); +} + +TEST_P (t8_test_cmesh_tree_reindex, small_bounding_box) +{ + check_lines ({ { 2e-10, 3e-10 }, { 0, 1e-10 }, { 1e-10, 2e-10 } }, { 2, 0, 1 }); +} + +TEST_P (t8_test_cmesh_tree_reindex, coincident_centers) +{ + /* Trees 0, 2, and 3 have the same center. Keep all of them in their original ID order. */ + check_lines ({ { 2, 4 }, { 0, 2 }, { 2.5, 3.5 }, { 1, 5 }, { 4, 6 } }, { 1, 0, 2, 3, 4 }); +} + +TEST_P (t8_test_cmesh_tree_reindex, single_tree) +{ + check_lines ({ { 0, 1 } }, { 0 }); +} + +TEST_P (t8_test_cmesh_tree_reindex, hypercube_preserves_geometry_and_connectivity) +{ + /* Periodic joins also exercise nonzero face orientations. */ + t8_cmesh_new_hypercube (&control_cmesh, T8_ECLASS_TET, sc_MPI_COMM_SELF, 0, 0, 1); + t8_cmesh_new_hypercube (&cmesh, T8_ECLASS_TET, sc_MPI_COMM_SELF, 0, 0, 1); + ASSERT_EQ (t8_cmesh_get_num_trees (cmesh), 6); + + /* Reconstruct the actual permutation independently of the retained map. */ + using vertex_array_t = std::array; + std::map vertices_to_old_id; + for (t8_locidx_t old_id = 0; old_id < 6; ++old_id) { + const double *vertices = t8_cmesh_get_tree_vertices (control_cmesh, old_id); + ASSERT_NE (vertices, nullptr); + vertex_array_t key; + std::copy_n (vertices, key.size (), key.begin ()); + ASSERT_TRUE (vertices_to_old_id.emplace (key, old_id).second); + } + std::vector observed_reindex (6, -1); + bool changed = false; + for (t8_locidx_t new_id = 0; new_id < 6; ++new_id) { + const double *vertices = t8_cmesh_get_tree_vertices (cmesh, new_id); + ASSERT_NE (vertices, nullptr); + vertex_array_t key; + std::copy_n (vertices, key.size (), key.begin ()); + const auto old_tree = vertices_to_old_id.find (key); + ASSERT_NE (old_tree, vertices_to_old_id.end ()); + observed_reindex[old_tree->second] = new_id; + changed = changed || old_tree->second != new_id; + } + EXPECT_EQ (changed, GetParam () >= 0); + check_mesh (observed_reindex); +} + +TEST_P (t8_test_cmesh_tree_reindex, partitioned_hypercube) +{ + check_partitioned_hypercube (T8_ECLASS_TET, 1); +} + +TEST_P (t8_test_cmesh_tree_reindex, partitioned_single_tree) +{ + check_partitioned_hypercube (T8_ECLASS_HEX, 0); +} + +TEST_P (t8_test_cmesh_tree_reindex, partitioned_distributed_lines) +{ + check_distributed_lines (false); +} + +TEST_P (t8_test_cmesh_tree_reindex, partitioned_lines_empty_root) +{ + check_distributed_lines (true); +} + +TEST_P (t8_test_cmesh_tree_reindex, partitioned_rank_local_map_retention) +{ + int mpirank; + int mpisize; + SC_CHECK_MPI (sc_MPI_Comm_rank (sc_MPI_COMM_WORLD, &mpirank)); + SC_CHECK_MPI (sc_MPI_Comm_size (sc_MPI_COMM_WORLD, &mpisize)); + const bool retain_map = GetParam () == 1 && mpirank == mpisize - 1; + if (GetParam () >= 0) { + t8_cmesh_enable_tree_reordering (cmesh, retain_map ? -1 : 0); + } + t8_cmesh_new_hypercube (&cmesh, T8_ECLASS_HEX, sc_MPI_COMM_WORLD, 0, 1, 0); + ASSERT_TRUE (t8_cmesh_is_committed (cmesh)); + EXPECT_EQ (t8_cmesh_get_num_trees (cmesh), 1); + EXPECT_EQ (t8_cmesh_get_num_local_trees (cmesh), mpirank == mpisize - 1 ? 1 : 0); + EXPECT_EQ (t8_cmesh_get_num_ghosts (cmesh), 0); + const t8_gloidx_t *reindex_map = t8_cmesh_get_tree_reindex_map (cmesh); + if (retain_map) { + ASSERT_NE (reindex_map, nullptr); + EXPECT_EQ (reindex_map[0], 0); + } + else { + EXPECT_EQ (reindex_map, nullptr); + } + if (mpirank == mpisize - 1) { + EXPECT_EQ (t8_cmesh_get_global_id (cmesh, 0), 0); + for (int iface = 0; iface < t8_eclass_num_faces[T8_ECLASS_HEX]; ++iface) { + EXPECT_TRUE (t8_cmesh_tree_face_is_boundary (cmesh, 0, iface)); + } + } + ASSERT_EQ (t8_cmesh_get_tree_reindex_map (cmesh), reindex_map); + if (reindex_map != nullptr) { + EXPECT_EQ (reindex_map[0], 0); + } +} + +INSTANTIATE_TEST_SUITE_P (reordering_modes, t8_test_cmesh_tree_reindex, testing::Values (-1, 0, 1), + [] (const testing::TestParamInfo &info) { + return info.param < 0 ? "disabled" : (info.param == 0 ? "discard_map" : "retain_map"); + });