diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8f7953ba1d..ea419d0375 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -151,6 +151,7 @@ target_sources( T8 PRIVATE t8_cmesh/t8_cmesh_internal/t8_cmesh_partition.cxx t8_cmesh/t8_cmesh_internal/t8_cmesh_stash.c 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 27c7e7cd79..a77ea8de5f 100644 --- a/src/t8_cmesh/t8_cmesh.cxx +++ b/src/t8_cmesh/t8_cmesh.cxx @@ -152,6 +152,12 @@ t8_cmesh_disable_negative_volume_check ([[maybe_unused]] t8_cmesh_t cmesh) #endif } +void +t8_cmesh_enable_tree_reordering (t8_cmesh_t cmesh) +{ + cmesh->reindex_trees = 1; +} + #if T8_ENABLE_DEBUG int t8_cmesh_validate_geometry (const t8_cmesh_t cmesh, const int check_for_negative_volume) @@ -239,6 +245,7 @@ t8_cmesh_init (t8_cmesh_t *pcmesh) #if T8_ENABLE_DEBUG cmesh->negative_volume_check = 1; #endif /* T8_ENABLE_DEBUG */ + cmesh->reindex_trees = 0; T8_ASSERT (t8_cmesh_is_initialized (cmesh)); } diff --git a/src/t8_cmesh/t8_cmesh.h b/src/t8_cmesh/t8_cmesh.h index 71c6520fe2..393711bc53 100644 --- a/src/t8_cmesh/t8_cmesh.h +++ b/src/t8_cmesh/t8_cmesh.h @@ -132,6 +132,13 @@ t8_cmesh_stash_is_empty (const t8_cmesh_t cmesh); void t8_cmesh_disable_negative_volume_check (t8_cmesh_t cmesh); +/** + * Enable localitly based indexing of trees during \ref t8_cmesh_commit. + * \param [in, out] cmesh + */ +void +t8_cmesh_enable_tree_reordering (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 +147,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 7d0f789302..9592c7e65e 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 /** * A struct to hold the information about a ghost facejoin. @@ -148,6 +149,10 @@ t8_cmesh_commit_replicated_new (t8_cmesh_t cmesh) t8_stash_class_struct_t *entry; t8_locidx_t num_trees = class_entries->elem_count, ltree; + if (cmesh->reindex_trees) { + t8_cmesh_tree_perform_reindex_inplace (stash, t8_cmesh_reindex_tree (cmesh)); + } + t8_cmesh_trees_init (&cmesh->trees, 1, num_trees, 0); t8_cmesh_trees_start_part (cmesh->trees, 0, 0, num_trees, 0, 0, 1); /* set tree classes */ 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..af47657d76 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -0,0 +1,399 @@ +/* + 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 + +/** \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. + */ + explicit t8_cmesh_tree_reindex_forest_data (const t8_locidx_t num_elements) + : elements (T8_ALLOC (t8_cmesh_tree_reindex_element_data, num_elements)), num_elements (num_elements), + 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. */ + 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. Unused. + * \param [in] lelement_id Tree-local leaf index of the considered element in \a forest_from. + * \param [in] scheme Refinement scheme of the forest. Unused. + * \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, + [[maybe_unused]] const t8_eclass_t tree_class, const t8_locidx_t lelement_id, + [[maybe_unused]] 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; + return data->elements[element_index].tree_centers.size () > 1 ? 1 : 0; +} + +/** + * \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. Unused. + * \param [in] scheme Refinement scheme of the forest. Unused. + * \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, + [[maybe_unused]] const t8_eclass_t tree_class, [[maybe_unused]] 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 (old_element_data.tree_centers.size () == old_element_data.tree_ids.size ()); + + if (refine == 0) { + /* An unchanged element can contain at most one center, otherwise the adapt callback would have refined it. */ + T8_ASSERT (num_old_elements == 1); + T8_ASSERT (num_new_elements == 1); + T8_ASSERT (old_element_data.tree_centers.size () <= 1); + + 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; + } + + T8_ASSERT (num_old_elements == 1); + 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 (), 0); + + /* 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; + } + } + + /* Any child that still contains multiple centers has to be refined in the next pass. */ + if (new_element_data.tree_centers.size () > 1) { + 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_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) +{ + 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 = { 0.0, 0.0, 0.0 }; + + 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); + } + + 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]; + constexpr double tolerance = T8_PRECISION_SQRT_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)); + + 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_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)); + 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 ()); + T8_ASSERT (element_data.tree_ids.size () <= 1); + + if (element_data.tree_ids.empty ()) { + continue; + } + + /* Non-empty leaves contain exactly one original tree ID at this point. */ + tree_reindex.emplace (element_data.tree_ids.front (), new_tree_index); + ++new_tree_index; + } + } + + 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); +} 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..90b6eedb1c --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.hxx @@ -0,0 +1,64 @@ +/* + 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. +*/ + +#ifndef T8_CMESH_REINDEX_TREES_H +#define T8_CMESH_REINDEX_TREES_H + +#include +#include +#include +#include + +/** + * Compute a spatially coherent reindexing of the coarse mesh trees. + * + * Constructs an auxiliary bounding box forest containing the centers of the + * original coarse mesh trees. The forest is adaptively refined until each + * leaf element contains at most one tree center. The leaves are then traversed + * in SFC order to assign new indices to the original trees. + * + * \param [in] cmesh The original coarse mesh whose trees are to be reindexed. + * \param [in] comm MPI communicator used to construct the auxiliary forest. + * + * \return A map from the original global tree IDs to their new SFC-based indices. + * + * \note This function only computes the reindexing and does not modify the + * original coarse mesh. + */ + +std::map +t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm = sc_MPI_COMM_SELF); + +/** + * Apply a tree reindexing to the coarse mesh stash in place. + * + * Updates the tree IDs stored in the stash according to the given mapping, + * allowing the coarse mesh to use the newly computed tree ordering. + * + * \param [in,out] stash The coarse mesh stash to be modified. + * \param [in] tree_reindex Mapping from original global tree IDs to + * their new indices. + */ +void +t8_cmesh_tree_perform_reindex_inplace (t8_stash_t &stash, const std::map &tree_reindex); + +#endif /* !T8_CMESH_REINDEX_TREES_H */ 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 b7e47b1c41..952377b615 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_types.h +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_types.h @@ -99,6 +99,8 @@ 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 */ + 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 dd8f4cc953..fe0556ef60 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -135,6 +135,7 @@ 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_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..b7844b53a1 --- /dev/null +++ b/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx @@ -0,0 +1,263 @@ +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +struct t8_test_cmesh_tree_reindex: public testing::Test +{ + protected: + static constexpr t8_gloidx_t num_trees = 6; + static constexpr int num_vertices = 4; + static constexpr int num_coords = 3 * num_vertices; + static constexpr int num_tet_faces = 4; + + using vertex_array_t = std::array; + + void + SetUp () override + { + t8_cmesh_init (&cmesh); + ASSERT_NE (cmesh, nullptr); + + /* + * The reindex_trees flag is initialized to 1 by default. + * Therefore, the test cmesh will be reindexed automatically during commit. + */ + if (!cmesh->reindex_trees) { + t8_cmesh_enable_tree_reordering (cmesh); + } + ASSERT_EQ (cmesh->reindex_trees, 1); + + t8_cmesh_init (&control_cmesh); + ASSERT_NE (control_cmesh, nullptr); + + /* + * The control mesh uses the same hypercube example, but disables reindexing. + * It is used as the reference mesh with the original tree order. + */ + control_cmesh->reindex_trees = 0; + } + + void + TearDown () override + { + if (cmesh != nullptr) { + t8_cmesh_unref (&cmesh); + cmesh = nullptr; + } + + if (control_cmesh != nullptr) { + t8_cmesh_unref (&control_cmesh); + control_cmesh = nullptr; + } + } + + static vertex_array_t + get_tree_vertices (t8_cmesh_t cmesh, const t8_locidx_t local_tree_id) + { + vertex_array_t vertices {}; + + double *tree_vertices = t8_cmesh_get_tree_vertices (cmesh, local_tree_id); + T8_ASSERT (tree_vertices != nullptr); + + for (int icoord = 0; icoord < num_coords; ++icoord) { + vertices[icoord] = tree_vertices[icoord]; + } + + return vertices; + } + + t8_cmesh_t cmesh = nullptr; + t8_cmesh_t control_cmesh = nullptr; +}; + +TEST_F (t8_test_cmesh_tree_reindex, hypercube_commit_reindexes_trees_successfully_and_correctly) +{ + t8_productionf ("Test started\n"); + + /* + * Construct the control cmesh from the existing t8code hypercube example. + * Reindexing was disabled in SetUp(), so this mesh keeps the original tree + * order generated by t8_cmesh_new_hypercube. + */ + t8_cmesh_new_hypercube (&control_cmesh, T8_ECLASS_TET, sc_MPI_COMM_SELF, 0, 0, 0); + + ASSERT_TRUE (t8_cmesh_is_committed (control_cmesh)); + + t8_cmesh_vtk_write_file (control_cmesh, "test_cmesh_tree_reindex_original"); + + /* + * Construct the actual test cmesh from the same hypercube example. + * The reindex_trees flag is enabled by default, so the reindexing is executed + * during the commit inside t8_cmesh_new_hypercube. + */ + t8_cmesh_new_hypercube (&cmesh, T8_ECLASS_TET, sc_MPI_COMM_SELF, 0, 0, 0); + + ASSERT_TRUE (t8_cmesh_is_committed (cmesh)); + + t8_cmesh_vtk_write_file (cmesh, "test_cmesh_tree_reindex_reindexed"); + + /* + * The tetrahedral hypercube consists of six tetrahedral trees. + * Since the test uses MPI_COMM_SELF, all trees are local. + */ + for (t8_gloidx_t tree_id = 0; tree_id < num_trees; ++tree_id) { + EXPECT_GE (t8_cmesh_get_local_id (control_cmesh, tree_id), 0); + EXPECT_GE (t8_cmesh_get_local_id (cmesh, tree_id), 0); + } + + /* + * Build a reference map from vertex coordinates to original tree ids. + * Since reindexing must not change the geometry, each reindexed tree should + * match exactly one tree of the control mesh. + */ + std::map vertices_to_old_tree_id; + + for (t8_gloidx_t old_tree_id = 0; old_tree_id < num_trees; ++old_tree_id) { + const t8_locidx_t old_local_tree_id = t8_cmesh_get_local_id (control_cmesh, old_tree_id); + + ASSERT_GE (old_local_tree_id, 0); + EXPECT_EQ (t8_cmesh_get_global_id (control_cmesh, old_local_tree_id), old_tree_id); + EXPECT_EQ (t8_cmesh_get_tree_class (control_cmesh, old_local_tree_id), T8_ECLASS_TET); + + const vertex_array_t vertices = get_tree_vertices (control_cmesh, old_local_tree_id); + + const auto inserted = vertices_to_old_tree_id.emplace (vertices, old_tree_id); + ASSERT_TRUE (inserted.second) << "Duplicate vertex data in control cmesh for tree " << old_tree_id; + } + + /* + * Reconstruct the observed reindexing map from the committed reindexed cmesh. + * If a tree in the reindexed cmesh has the same vertex coordinates as an old + * tree in the control cmesh, then this old tree was mapped to the current new + * tree id. + */ + std::map observed_reindex; + std::set old_tree_ids; + std::set new_tree_ids; + + bool reindex_is_identity = true; + + for (t8_gloidx_t new_tree_id = 0; new_tree_id < num_trees; ++new_tree_id) { + const t8_locidx_t new_local_tree_id = t8_cmesh_get_local_id (cmesh, new_tree_id); + + ASSERT_GE (new_local_tree_id, 0); + EXPECT_EQ (t8_cmesh_get_global_id (cmesh, new_local_tree_id), new_tree_id); + EXPECT_EQ (t8_cmesh_get_tree_class (cmesh, new_local_tree_id), T8_ECLASS_TET); + + const vertex_array_t vertices = get_tree_vertices (cmesh, new_local_tree_id); + + const auto old_tree_entry = vertices_to_old_tree_id.find (vertices); + + ASSERT_NE (old_tree_entry, vertices_to_old_tree_id.end ()) + << "Could not find matching original tree for reindexed tree " << new_tree_id; + + const t8_gloidx_t old_tree_id = old_tree_entry->second; + + const auto inserted = observed_reindex.emplace (old_tree_id, new_tree_id); + ASSERT_TRUE (inserted.second) << "Old tree " << old_tree_id << " was mapped more than once."; + + old_tree_ids.insert (old_tree_id); + new_tree_ids.insert (new_tree_id); + + if (old_tree_id != new_tree_id) { + reindex_is_identity = false; + } + + t8_productionf ("Observed reindex: old global tree id %lli -> new global tree id %lli\n", + static_cast (old_tree_id), static_cast (new_tree_id)); + } + + /* + * Check that the observed reindexing is a valid bijection. + */ + EXPECT_EQ (observed_reindex.size (), static_cast (num_trees)); + EXPECT_EQ (old_tree_ids.size (), static_cast (num_trees)); + EXPECT_EQ (new_tree_ids.size (), static_cast (num_trees)); + + for (t8_gloidx_t tree_id = 0; tree_id < num_trees; ++tree_id) { + EXPECT_EQ (old_tree_ids.count (tree_id), 1); + EXPECT_EQ (new_tree_ids.count (tree_id), 1); + } + + t8_productionf ("Reindexing is %s\n", reindex_is_identity ? "identity" : "non-identity"); + + /* + * Verify that all internal face joins were updated consistently. + * The control cmesh gives the original connectivity. For each original join, + * the old tree ids are mapped to their new tree ids and the same face + * connection is queried in the reindexed cmesh. + */ + for (t8_gloidx_t old_tree_id = 0; old_tree_id < num_trees; ++old_tree_id) { + const t8_locidx_t old_local_tree_id = t8_cmesh_get_local_id (control_cmesh, old_tree_id); + + ASSERT_GE (old_local_tree_id, 0); + + for (int face = 0; face < num_tet_faces; ++face) { + int old_dual_face = -1; + int old_orientation = -1; + + const t8_locidx_t old_neighbor_local_tree_id + = t8_cmesh_get_face_neighbor (control_cmesh, old_local_tree_id, face, &old_dual_face, &old_orientation); + + if (old_neighbor_local_tree_id < 0) { + continue; + } + + const t8_gloidx_t old_neighbor_tree_id = t8_cmesh_get_global_id (control_cmesh, old_neighbor_local_tree_id); + + /* + * Check every undirected join only once. + */ + if (old_tree_id > old_neighbor_tree_id) { + continue; + } + + const t8_gloidx_t new_tree_id = observed_reindex.at (old_tree_id); + const t8_gloidx_t new_neighbor_tree_id = observed_reindex.at (old_neighbor_tree_id); + + const t8_locidx_t new_local_tree_id = t8_cmesh_get_local_id (cmesh, new_tree_id); + const t8_locidx_t new_neighbor_local_tree_id = t8_cmesh_get_local_id (cmesh, new_neighbor_tree_id); + + ASSERT_GE (new_local_tree_id, 0); + ASSERT_GE (new_neighbor_local_tree_id, 0); + + int new_dual_face = -1; + int new_orientation = -1; + + const t8_locidx_t actual_neighbor_local_tree_id + = t8_cmesh_get_face_neighbor (cmesh, new_local_tree_id, face, &new_dual_face, &new_orientation); + + ASSERT_GE (actual_neighbor_local_tree_id, 0); + + EXPECT_EQ (t8_cmesh_get_global_id (cmesh, actual_neighbor_local_tree_id), new_neighbor_tree_id); + EXPECT_EQ (new_dual_face, old_dual_face); + EXPECT_EQ (new_orientation, old_orientation); + + /* + * Also verify the reverse direction of the join. + */ + int reverse_dual_face = -1; + int reverse_orientation = -1; + + const t8_locidx_t reverse_neighbor_local_tree_id = t8_cmesh_get_face_neighbor ( + cmesh, new_neighbor_local_tree_id, old_dual_face, &reverse_dual_face, &reverse_orientation); + + ASSERT_GE (reverse_neighbor_local_tree_id, 0); + + EXPECT_EQ (t8_cmesh_get_global_id (cmesh, reverse_neighbor_local_tree_id), new_tree_id); + EXPECT_EQ (reverse_dual_face, face); + + t8_productionf ("Verified reindexed join old trees %lli-%lli -> new trees %lli-%lli\n", + static_cast (old_tree_id), static_cast (old_neighbor_tree_id), + static_cast (new_tree_id), static_cast (new_neighbor_tree_id)); + } + } +}