From fdea0bfd169facd42a0c90aa7cd75ef4e0a23f7b Mon Sep 17 00:00:00 2001 From: tunxu Date: Fri, 15 May 2026 17:11:46 +0200 Subject: [PATCH 01/10] initial Commit --- .../t8_cmesh_tree_reindex.cxx | 26 +++++++++++++++++ .../t8_cmesh_tree_reindex.hxx | 29 +++++++++++++++++++ 2 files changed, 55 insertions(+) create mode 100644 src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx create mode 100644 src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.hxx 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..1a7df86cb1 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -0,0 +1,26 @@ +/* + 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. +*/ +#include + +std::map t8_cmesh_reindex_tree(t8_cmesh_t cmesh){ + +}; \ No newline at end of file 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..230c85b238 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.hxx @@ -0,0 +1,29 @@ +/* + 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. +*/ + +#include +#include +#include +#include + +std::map t8_cmesh_reindex_tree(t8_cmesh_t cmesh){}; + From 0ab2e075e75ac3e89c8e5ded29219770db534468 Mon Sep 17 00:00:00 2001 From: tunxu Date: Tue, 19 May 2026 14:51:26 +0200 Subject: [PATCH 02/10] Create a cmesh from bounding box around given cmesh --- .../t8_cmesh_tree_reindex.cxx | 46 +++++++++++++++++-- 1 file changed, 43 insertions(+), 3 deletions(-) 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 index 1a7df86cb1..8159c28326 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -20,7 +20,47 @@ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ #include +#include +#include +#include +#include -std::map t8_cmesh_reindex_tree(t8_cmesh_t cmesh){ - -}; \ No newline at end of file +std::array +bbox_bounds_to_hex_vertices (const double bounds[6]) +{ + const double x_min = bounds[0]; + const double x_max = bounds[1]; + const double y_min = bounds[2]; + const double y_max = bounds[3]; + const double z_min = bounds[4]; + const double z_max = bounds[5]; + + return { + x_min, y_min, z_min, // vertex 0 + x_max, y_min, z_min, // vertex 1 + x_min, y_max, z_min, // vertex 2 + x_max, y_max, z_min, // vertex 3 + + x_min, y_min, z_max, // vertex 4 + x_max, y_min, z_max, // vertex 5 + x_min, y_max, z_max, // vertex 6 + x_max, y_max, z_max // vertex 7 + }; +} + +std::map +t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) +{ + double bounding_box[6]; + t8_cmesh_get_local_bounding_box (cmesh, bounding_box); + + t8_cmesh_t bbox_cmesh; + + t8_cmesh_init (&bbox_cmesh); + t8_cmesh_set_tree_class (bbox_cmesh, 0, T8_ECLASS_HEX); + std::array vertices = bbox_bounds_to_hex_vertices (bounding_box); + + t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 8); + + t8_cmesh_commit (bbox_cmesh, comm); +}; From 2c69750018d2e31149df3d7fe5cb5623b8a2b11a Mon Sep 17 00:00:00 2001 From: tunxu Date: Tue, 26 May 2026 10:56:50 +0200 Subject: [PATCH 03/10] Added check for required refinement level --- src/CMakeLists.txt | 1 + .../t8_cmesh_tree_reindex.cxx | 240 +++++++++++++++++- .../t8_cmesh_tree_reindex.hxx | 4 +- 3 files changed, 236 insertions(+), 9 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 53a710f213..d81d7c3e25 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -152,6 +152,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_internal/t8_cmesh_tree_reindex.cxx b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx index 8159c28326..6e3b80fb63 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -19,11 +19,20 @@ along with t8code; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ -#include + +#include + #include #include #include +#include + #include +#include +#include +#include +#include +#include std::array bbox_bounds_to_hex_vertices (const double bounds[6]) @@ -48,19 +57,236 @@ bbox_bounds_to_hex_vertices (const double bounds[6]) }; } -std::map +std::array +get_reference_center (t8_eclass_t eclass) +{ + switch (eclass) { + case T8_ECLASS_TRIANGLE: + return { 1.0 / 3.0, 1.0 / 3.0, 0.0 }; + + case T8_ECLASS_QUAD: + return { 0.5, 0.5, 0.0 }; + + case T8_ECLASS_HEX: + return { 0.5, 0.5, 0.5 }; + + case T8_ECLASS_TET: + return { 0.25, 0.25, 0.25 }; + + default: + return { 0.5, 0.5, 0.5 }; + } +} + +std::array +compute_tree_midpoint_with_geometry (t8_forest_t forest, t8_locidx_t local_tree_id) +{ + const t8_eclass_t eclass = t8_forest_get_eclass (forest, local_tree_id); + + const std::array ref_midpoint = get_reference_center (eclass); + + double physical_midpoint[3] = { 0.0, 0.0, 0.0 }; + + const t8_element_t *element = t8_forest_get_leaf_element (forest, local_tree_id, 0); + + t8_forest_element_from_ref_coords (forest, local_tree_id, element, ref_midpoint.data (), 1, physical_midpoint); + + return { physical_midpoint[0], physical_midpoint[1], physical_midpoint[2] }; +} + +struct cell_coordinate +{ + t8_locidx_t x; + t8_locidx_t y; + t8_locidx_t z; +}; + +cell_coordinate +point_to_bbox_cell (const std::array &point, const double bounds[6], const int level) +{ + const t8_locidx_t cells_per_direction = static_cast (1) << level; + + const double x_min = bounds[0]; + const double x_max = bounds[1]; + const double y_min = bounds[2]; + const double y_max = bounds[3]; + const double z_min = bounds[4]; + const double z_max = bounds[5]; + + auto compute_index + = [cells_per_direction] (const double value, const double min_value, const double max_value) -> t8_locidx_t { + const double length = max_value - min_value; + + if (length <= 0.0) { + return 0; + } + + const double normalized = (value - min_value) / length; + + t8_locidx_t index = static_cast (std::floor (normalized * cells_per_direction)); + + if (index < 0) { + index = 0; + } + + if (index >= cells_per_direction) { + index = cells_per_direction - 1; + } + + return index; + }; + + return { compute_index (point[0], x_min, x_max), compute_index (point[1], y_min, y_max), + compute_index (point[2], z_min, z_max) }; +} + +bool +refinement_level_is_unique (const std::map> &tree_to_center, const double bounds[6], + const int level) +{ + std::set> occupied_cells; + + for (const auto &entry : tree_to_center) { + const std::array ¢er = entry.second; + + const cell_coordinate cell = point_to_bbox_cell (center, bounds, level); + + const auto key = std::make_tuple (cell.x, cell.y, cell.z); + + const auto insert_result = occupied_cells.insert (key); + + if (!insert_result.second) { + return false; + } + } + + return true; +} + +int +find_required_refinement_level (const std::map> &tree_to_center, + const double bounds[6], const int max_level) +{ + for (int level = 0; level <= max_level; ++level) { + if (refinement_level_is_unique (tree_to_center, bounds, level)) { + return level; + } + } + + return -1; +} + +std::map t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) { + std::map tree_reindex; + double bounding_box[6]; t8_cmesh_get_local_bounding_box (cmesh, bounding_box); - t8_cmesh_t bbox_cmesh; + auto vertices = bbox_bounds_to_hex_vertices (bounding_box); + /* + * Build auxiliary bounding-box cmesh. + */ + t8_cmesh_t bbox_cmesh; t8_cmesh_init (&bbox_cmesh); - t8_cmesh_set_tree_class (bbox_cmesh, 0, T8_ECLASS_HEX); - std::array vertices = bbox_bounds_to_hex_vertices (bounding_box); + t8_cmesh_set_tree_class (bbox_cmesh, 0, T8_ECLASS_HEX); t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 8); - t8_cmesh_commit (bbox_cmesh, comm); -}; + sc_MPI_Comm bbox_comm = comm; + + t8_cmesh_commit (bbox_cmesh, bbox_comm); + + /* + * Build a level-0 forest from the original cmesh. + * This is only used to evaluate the physical centers of the original trees. + * + * t8_forest_new_uniform takes ownership of a cmesh reference. + * Since cmesh belongs to the caller, increase the refcount first. + */ + t8_cmesh_ref (cmesh); + + t8_forest_t center_forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default (), 0, 0, comm); + + const t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (center_forest); + + /* + * Compute one center point per original local tree. + */ + std::map> tree_to_center; + + for (t8_locidx_t itree = 0; itree < num_local_trees; ++itree) { + std::array midpoint = compute_tree_midpoint_with_geometry (center_forest, itree); + + tree_to_center.insert ({ itree, midpoint }); + + t8_productionf ("Tree Index: %u \t Tree Center: %f, %f, %f\n", static_cast (itree), midpoint[0], + midpoint[1], midpoint[2]); + } + + /* + * Now find the first bounding-box refinement level where each cell contains + * at most one tree center. + */ + const int max_level = 29; + + const int required_level = find_required_refinement_level (tree_to_center, bounding_box, max_level); + + if (required_level < 0) { + t8_productionf ("Could not find a refinement level where all tree centers are separated.\n" + "This can happen if two trees have exactly the same center.\n"); + + /* + * Fallback: identity reindexing. + */ + for (t8_locidx_t itree = 0; itree < num_local_trees; ++itree) { + tree_reindex.insert ({ itree, itree }); + } + + t8_forest_unref (¢er_forest); + + /* + * bbox_cmesh was committed but not consumed by a forest in this branch. + */ + t8_cmesh_unref (&bbox_cmesh); + + return tree_reindex; + } + + t8_productionf ("Required bounding-box refinement level: %i\n", required_level); + + /* + * Now create the actual refined bounding-box forest. + * This forest represents your refined auxiliary bounding box. + * + * t8_forest_new_uniform takes ownership of bbox_cmesh and of the scheme. + */ + t8_forest_t bbox_forest = t8_forest_new_uniform (bbox_cmesh, t8_scheme_new_default (), required_level, 0, bbox_comm); + + /* + * At this point: + * + * - center_forest stores the original cmesh trees at level 0. + * - bbox_forest stores the refined bounding box. + * - required_level satisfies your criterion. + * + * For now, fill identity reindexing. + * The next step will be to replace this with SFC/Morton ordering. + */ + for (t8_locidx_t itree = 0; itree < num_local_trees; ++itree) { + const cell_coordinate cell = point_to_bbox_cell (tree_to_center[itree], bounding_box, required_level); + + tree_reindex.insert ({ itree, itree }); + + t8_productionf ("Tree %u lies in bbox cell (%u, %u, %u) at level %i\n", static_cast (itree), + static_cast (cell.x), static_cast (cell.y), static_cast (cell.z), + required_level); + } + + t8_forest_unref (&bbox_forest); + t8_forest_unref (¢er_forest); + + return tree_reindex; +} 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 index 230c85b238..89d345e016 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.hxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.hxx @@ -25,5 +25,5 @@ #include #include -std::map t8_cmesh_reindex_tree(t8_cmesh_t cmesh){}; - +std::map +t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm); From a480c714d35e6fc9dbba263376f3a3ab746e0fdd Mon Sep 17 00:00:00 2001 From: tunxu Date: Wed, 3 Jun 2026 17:42:33 +0200 Subject: [PATCH 04/10] working test --- .../t8_cmesh_tree_reindex.cxx | 155 ++++++++++++++---- test/CMakeLists.txt | 1 + test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx | 74 +++++++++ 3 files changed, 197 insertions(+), 33 deletions(-) create mode 100644 test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx 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 index 6e3b80fb63..23d5c029e3 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -20,19 +20,26 @@ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ +/* + 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. +*/ + #include #include #include #include #include +#include #include -#include -#include #include +#include #include #include +#include std::array bbox_bounds_to_hex_vertices (const double bounds[6]) @@ -101,6 +108,12 @@ struct cell_coordinate t8_locidx_t z; }; +static std::tuple +cell_to_key (const cell_coordinate &cell) +{ + return std::make_tuple (cell.x, cell.y, cell.z); +} + cell_coordinate point_to_bbox_cell (const std::array &point, const double bounds[6], const int level) { @@ -151,7 +164,7 @@ refinement_level_is_unique (const std::map> & const cell_coordinate cell = point_to_bbox_cell (center, bounds, level); - const auto key = std::make_tuple (cell.x, cell.y, cell.z); + const auto key = cell_to_key (cell); const auto insert_result = occupied_cells.insert (key); @@ -176,35 +189,51 @@ find_required_refinement_level (const std::map (1) << level; + + return static_cast (cell.x) + cells_per_direction * static_cast (cell.y) + + cells_per_direction * cells_per_direction * static_cast (cell.z); +} + std::map t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) { std::map tree_reindex; + /* + * 1. Get local bounding box of the original cmesh. + */ double bounding_box[6]; t8_cmesh_get_local_bounding_box (cmesh, bounding_box); - auto vertices = bbox_bounds_to_hex_vertices (bounding_box); - /* - * Build auxiliary bounding-box cmesh. + * 2. Build auxiliary bounding-box cmesh. */ + const auto vertices = bbox_bounds_to_hex_vertices (bounding_box); + t8_cmesh_t bbox_cmesh; t8_cmesh_init (&bbox_cmesh); t8_cmesh_set_tree_class (bbox_cmesh, 0, T8_ECLASS_HEX); t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 8); + /* + * If bounding_box is process-local, SC_MPI_COMM_SELF is usually safer. + * If you intentionally want this bbox cmesh on the same communicator, + * keep bbox_comm = comm. + */ sc_MPI_Comm bbox_comm = comm; t8_cmesh_commit (bbox_cmesh, bbox_comm); /* - * Build a level-0 forest from the original cmesh. - * This is only used to evaluate the physical centers of the original trees. + * 3. Build level-0 forest from original cmesh to evaluate tree centers. * * t8_forest_new_uniform takes ownership of a cmesh reference. - * Since cmesh belongs to the caller, increase the refcount first. + * Since cmesh belongs to the caller, we increase the refcount first. */ t8_cmesh_ref (cmesh); @@ -213,12 +242,12 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) const t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (center_forest); /* - * Compute one center point per original local tree. + * 4. Compute one physical center point per original local tree. */ std::map> tree_to_center; for (t8_locidx_t itree = 0; itree < num_local_trees; ++itree) { - std::array midpoint = compute_tree_midpoint_with_geometry (center_forest, itree); + const std::array midpoint = compute_tree_midpoint_with_geometry (center_forest, itree); tree_to_center.insert ({ itree, midpoint }); @@ -227,8 +256,8 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) } /* - * Now find the first bounding-box refinement level where each cell contains - * at most one tree center. + * 5. Find the first bbox refinement level where each cell contains + * at most one tree center. */ const int max_level = 29; @@ -238,11 +267,8 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) t8_productionf ("Could not find a refinement level where all tree centers are separated.\n" "This can happen if two trees have exactly the same center.\n"); - /* - * Fallback: identity reindexing. - */ for (t8_locidx_t itree = 0; itree < num_local_trees; ++itree) { - tree_reindex.insert ({ itree, itree }); + tree_reindex[itree] = itree; } t8_forest_unref (¢er_forest); @@ -258,33 +284,96 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) t8_productionf ("Required bounding-box refinement level: %i\n", required_level); /* - * Now create the actual refined bounding-box forest. - * This forest represents your refined auxiliary bounding box. + * 6. Map each original tree center to the logical bbox-cell linear id + * at required_level. * - * t8_forest_new_uniform takes ownership of bbox_cmesh and of the scheme. + * Since required_level passed the uniqueness test, each logical cell + * should contain at most one tree center. */ - t8_forest_t bbox_forest = t8_forest_new_uniform (bbox_cmesh, t8_scheme_new_default (), required_level, 0, bbox_comm); + std::map linear_id_to_tree; + + for (const auto &entry : tree_to_center) { + const t8_locidx_t original_tree = entry.first; + const std::array ¢er = entry.second; + + const cell_coordinate cell = point_to_bbox_cell (center, bounding_box, required_level); + + const t8_linearidx_t linear_id = cell_to_linear_id (cell, required_level); + + const auto insert_result = linear_id_to_tree.insert ({ linear_id, original_tree }); + + if (!insert_result.second) { + /* + * This should not happen because required_level was checked. + */ + t8_productionf ("Unexpected collision: more than one center in cell (%u, %u, %u), linear id %li.\n", + static_cast (cell.x), static_cast (cell.y), static_cast (cell.z), + static_cast (linear_id)); + } + + t8_productionf ("Tree %u lies in bbox cell (%u, %u, %u), linear id %li at level %i\n", + static_cast (original_tree), static_cast (cell.x), + static_cast (cell.y), static_cast (cell.z), static_cast (linear_id), + required_level); + } /* - * At this point: + * 7. Create the refined bbox forest. * - * - center_forest stores the original cmesh trees at level 0. - * - bbox_forest stores the refined bounding box. - * - required_level satisfies your criterion. + * This forest is only used to get the t8code/SFC order of the bbox leaves. + * t8_forest_new_uniform takes ownership of bbox_cmesh and the scheme. + */ + t8_forest_t bbox_forest = t8_forest_new_uniform (bbox_cmesh, t8_scheme_new_default (), required_level, 0, bbox_comm); + + /* + * 8. Iterate over bbox leaves in t8code order. * - * For now, fill identity reindexing. - * The next step will be to replace this with SFC/Morton ordering. + * We do not evaluate geometry on bbox_forest here. + * Instead, we use each leaf's linear id at required_level. */ - for (t8_locidx_t itree = 0; itree < num_local_trees; ++itree) { - const cell_coordinate cell = point_to_bbox_cell (tree_to_center[itree], bounding_box, required_level); + t8_locidx_t next_sfc_index = 0; - tree_reindex.insert ({ itree, itree }); + const t8_locidx_t num_bbox_local_trees = t8_forest_get_num_local_trees (bbox_forest); - t8_productionf ("Tree %u lies in bbox cell (%u, %u, %u) at level %i\n", static_cast (itree), - static_cast (cell.x), static_cast (cell.y), static_cast (cell.z), - required_level); + for (t8_locidx_t bbox_tree = 0; bbox_tree < num_bbox_local_trees; ++bbox_tree) { + t8_element_array_t *leaf_elements = t8_forest_tree_get_leaf_elements (bbox_forest, bbox_tree); + + for (auto leaf_it = t8_element_array_begin (leaf_elements); leaf_it != t8_element_array_end (leaf_elements); + ++leaf_it) { + const t8_linearidx_t leaf_linear_id = leaf_it.get_linear_id_at_level (required_level); + + const auto found_tree = linear_id_to_tree.find (leaf_linear_id); + + if (found_tree == linear_id_to_tree.end ()) { + /* + * Empty bbox leaf. No original tree center lies here. + */ + continue; + } + + const t8_locidx_t original_tree = found_tree->second; + + tree_reindex[original_tree] = next_sfc_index; + + t8_productionf ("bbox leaf linear id %li maps original tree %u to new SFC index %u\n", + static_cast (leaf_linear_id), static_cast (original_tree), + static_cast (next_sfc_index)); + + ++next_sfc_index; + } } + /* + * 9. Sanity check. + */ + if (tree_reindex.size () != static_cast (num_local_trees)) { + t8_productionf ("Warning: only mapped %u of %u local trees.\n", static_cast (tree_reindex.size ()), + static_cast (num_local_trees)); + } + + /* + * 10. Cleanup. + */ t8_forest_unref (&bbox_forest); t8_forest_unref (¢er_forest); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 80ca77276d..e6b845d68a 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_compute_first_element_serial SOURCES t8_cmesh/t8_gtest_compute_first_element.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 ) 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..d68b7d520f --- /dev/null +++ b/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx @@ -0,0 +1,74 @@ +#include +#include +#include + +#include + +#include +#include +#include + +struct t8_test_cmesh_tree_reindex: public testing::Test +{ + protected: + void + SetUp () override + { + const int do_partition = 0; + const int periodic = 0; + + cmesh = t8_cmesh_new_hypercube (T8_ECLASS_TRIANGLE, sc_MPI_COMM_SELF, do_partition, periodic, 0); + + ASSERT_NE (cmesh, nullptr); + } + + void + TearDown () override + { + if (cmesh != nullptr) { + t8_cmesh_unref (&cmesh); + cmesh = nullptr; + } + } + + t8_cmesh_t cmesh = nullptr; +}; + +TEST_F (t8_test_cmesh_tree_reindex, sanity_check) +{ + t8_productionf ("Test started\n"); + + const t8_locidx_t num_local_trees = t8_cmesh_get_num_local_trees (cmesh); + + EXPECT_EQ (num_local_trees, 2); + + std::map reindex = t8_cmesh_reindex_tree (cmesh, sc_MPI_COMM_SELF); + + EXPECT_EQ (reindex.size (), static_cast (num_local_trees)); + + std::set new_indices; + + for (const auto &entry : reindex) { + const t8_locidx_t old_tree_id = entry.first; + const t8_locidx_t new_sfc_index = entry.second; + + EXPECT_GE (old_tree_id, 0); + EXPECT_LT (old_tree_id, num_local_trees); + + EXPECT_GE (new_sfc_index, 0); + EXPECT_LT (new_sfc_index, num_local_trees); + + new_indices.insert (new_sfc_index); + + t8_productionf ("old local tree id %li -> new local SFC index %li\n", static_cast (old_tree_id), + static_cast (new_sfc_index)); + } + + EXPECT_EQ (new_indices.size (), static_cast (num_local_trees)); + + /* + * Optional: for two trees, the SFC indices should be exactly {0, 1}. + */ + EXPECT_TRUE (new_indices.count (0) == 1); + EXPECT_TRUE (new_indices.count (1) == 1); +} From 06be9b5a4bfefa119830ef925055fea6771f4a90 Mon Sep 17 00:00:00 2001 From: tunxu Date: Mon, 15 Jun 2026 09:41:18 +0200 Subject: [PATCH 05/10] Zwischencommit zum Umwechseln --- .../t8_cmesh_internal/t8_cmesh_tree_reindex.cxx | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) 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 index 23d5c029e3..d994f04c21 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -220,11 +220,6 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) t8_cmesh_set_tree_class (bbox_cmesh, 0, T8_ECLASS_HEX); t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 8); - /* - * If bounding_box is process-local, SC_MPI_COMM_SELF is usually safer. - * If you intentionally want this bbox cmesh on the same communicator, - * keep bbox_comm = comm. - */ sc_MPI_Comm bbox_comm = comm; t8_cmesh_commit (bbox_cmesh, bbox_comm); @@ -320,7 +315,7 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) /* * 7. Create the refined bbox forest. * - * This forest is only used to get the t8code/SFC order of the bbox leaves. + * This forest is only used to get the SFC order of the bbox leaves. * t8_forest_new_uniform takes ownership of bbox_cmesh and the scheme. */ t8_forest_t bbox_forest = t8_forest_new_uniform (bbox_cmesh, t8_scheme_new_default (), required_level, 0, bbox_comm); @@ -379,3 +374,11 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) return tree_reindex; } + +//void +//t8_perform_reindex (t8_cmesh_t cmesh, std::map tree_reindex_map) +//{ +// int num_trees = cmesh->num_trees; +// for (t8_locidx_t itree = 0; itree < num_trees; itree++) { +// } +//} From 82e1524c4f57894b7b6adeb35082a4db864a272c Mon Sep 17 00:00:00 2001 From: Tu Nguyen Xuan Date: Wed, 24 Jun 2026 13:57:31 +0200 Subject: [PATCH 06/10] Fixed some issues --- .../t8_cmesh_tree_reindex.cxx | 553 +++++++++++------- 1 file changed, 333 insertions(+), 220 deletions(-) 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 index d994f04c21..876db975ef 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -12,190 +12,214 @@ 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. + 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. */ -/* - 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. -*/ - #include +#include +#include #include #include #include #include #include +#include +#include +#include +#include #include #include #include -#include -#include #include +#include -std::array -bbox_bounds_to_hex_vertices (const double bounds[6]) -{ - const double x_min = bounds[0]; - const double x_max = bounds[1]; - const double y_min = bounds[2]; - const double y_max = bounds[3]; - const double z_min = bounds[4]; - const double z_max = bounds[5]; - - return { - x_min, y_min, z_min, // vertex 0 - x_max, y_min, z_min, // vertex 1 - x_min, y_max, z_min, // vertex 2 - x_max, y_max, z_min, // vertex 3 - - x_min, y_min, z_max, // vertex 4 - x_max, y_min, z_max, // vertex 5 - x_min, y_max, z_max, // vertex 6 - x_max, y_max, z_max // vertex 7 - }; -} - -std::array -get_reference_center (t8_eclass_t eclass) -{ - switch (eclass) { - case T8_ECLASS_TRIANGLE: - return { 1.0 / 3.0, 1.0 / 3.0, 0.0 }; - - case T8_ECLASS_QUAD: - return { 0.5, 0.5, 0.0 }; - - case T8_ECLASS_HEX: - return { 0.5, 0.5, 0.5 }; +//static int +//bbox_is_effectively_2d (const double bounds[6]) +//{ +// const double dz = std::fabs (bounds[5] - bounds[4]); +// const double scale = std::max ({ 1.0, std::fabs (bounds[4]), std::fabs (bounds[5]) }); +// return dz <= 1e-14 * scale; +//} - case T8_ECLASS_TET: - return { 0.25, 0.25, 0.25 }; +//static double +//normalize_coordinate (const double value, const double min_value, const double max_value, const double tolerance) +//{ +// const double length = max_value - min_value; +// T8_ASSERT (std::fabs (length) > tolerance); +// +// double normalized = (value - min_value) / length; +// +// /* Avoid losing points due to tiny roundoff outside the unit bbox. */ +// if (normalized < 0.0 && normalized >= -tolerance) { +// normalized = 0.0; +// } +// if (normalized > 1.0 && normalized <= 1.0 + tolerance) { +// normalized = 1.0; +// } +// +// return normalized; +//} - default: - return { 0.5, 0.5, 0.5 }; - } -} +//static std::array +//unit_quad_vertices () +//{ +// return { +// 0.0, 0.0, 0.0, // vertex 0 +// 1.0, 0.0, 0.0, // vertex 1 +// 0.0, 1.0, 0.0, // vertex 2 +// 1.0, 1.0, 0.0 // vertex 3 +// }; +//} +// +//static std::array +//unit_hex_vertices () +//{ +// return { +// 0.0, 0.0, 0.0, // vertex 0 +// 1.0, 0.0, 0.0, // vertex 1 +// 0.0, 1.0, 0.0, // vertex 2 +// 1.0, 1.0, 0.0, // vertex 3 +// +// 0.0, 0.0, 1.0, // vertex 4 +// 1.0, 0.0, 1.0, // vertex 5 +// 0.0, 1.0, 1.0, // vertex 6 +// 1.0, 1.0, 1.0 // vertex 7 +// }; +//} std::array compute_tree_midpoint_with_geometry (t8_forest_t forest, t8_locidx_t local_tree_id) { - const t8_eclass_t eclass = t8_forest_get_eclass (forest, local_tree_id); - - const std::array ref_midpoint = get_reference_center (eclass); - - double physical_midpoint[3] = { 0.0, 0.0, 0.0 }; + std::array physical_midpoint = { 0.0, 0.0, 0.0 }; const t8_element_t *element = t8_forest_get_leaf_element (forest, local_tree_id, 0); - t8_forest_element_from_ref_coords (forest, local_tree_id, element, ref_midpoint.data (), 1, physical_midpoint); + t8_forest_element_centroid (forest, local_tree_id, element, physical_midpoint.data ()); - return { physical_midpoint[0], physical_midpoint[1], physical_midpoint[2] }; -} - -struct cell_coordinate -{ - t8_locidx_t x; - t8_locidx_t y; - t8_locidx_t z; -}; + t8_productionf ("Computed centroid for local tree %u: %.17g, %.17g, %.17g\n", static_cast (local_tree_id), + physical_midpoint[0], physical_midpoint[1], physical_midpoint[2]); -static std::tuple -cell_to_key (const cell_coordinate &cell) -{ - return std::make_tuple (cell.x, cell.y, cell.z); + return physical_midpoint; } -cell_coordinate -point_to_bbox_cell (const std::array &point, const double bounds[6], const int level) +/* + * Data passed to the adapt callback. + * + * points_flat stores all normalized original tree centers as: + * + * x0, y0, z0, x1, y1, z1, ... + * + * The helper bbox cmesh is deliberately built as a unit quad/hex. Therefore + * these points are normalized into [0, 1]^d before they are passed to + * t8_forest_element_points_inside. + */ +struct bbox_adapt_data { - const t8_locidx_t cells_per_direction = static_cast (1) << level; - - const double x_min = bounds[0]; - const double x_max = bounds[1]; - const double y_min = bounds[2]; - const double y_max = bounds[3]; - const double z_min = bounds[4]; - const double z_max = bounds[5]; - - auto compute_index - = [cells_per_direction] (const double value, const double min_value, const double max_value) -> t8_locidx_t { - const double length = max_value - min_value; - - if (length <= 0.0) { - return 0; - } - - const double normalized = (value - min_value) / length; - - t8_locidx_t index = static_cast (std::floor (normalized * cells_per_direction)); - - if (index < 0) { - index = 0; - } + const double *points_flat; + int num_points; + double tolerance; + int max_level; - if (index >= cells_per_direction) { - index = cells_per_direction - 1; - } - - return index; - }; - - return { compute_index (point[0], x_min, x_max), compute_index (point[1], y_min, y_max), - compute_index (point[2], z_min, z_max) }; -} + /* Set to 1 by the callback whenever at least one leaf is refined. */ + //int refined_any; +}; -bool -refinement_level_is_unique (const std::map> &tree_to_center, const double bounds[6], - const int level) +/* + * Count how many original tree centers are inside a given bbox leaf. + * + * This intentionally uses the t8code implementation t8_forest_element_points_inside. + * The points must be in the coordinate system of the helper bbox forest. Since + * the helper bbox forest is the unit quad/hex, points_flat contains normalized + * coordinates, not physical coordinates. + * + * If contained_tree_id is not nullptr and exactly one point is inside, + * contained_tree_id will be set to the corresponding original tree id. + */ +static int +count_tree_centers_inside_leaf1 (t8_forest_t forest, t8_locidx_t which_tree, const t8_element_t *element, + const bbox_adapt_data *data, t8_locidx_t *contained_tree_id) { - std::set> occupied_cells; + std::vector points_inside (static_cast (data->num_points), 0); - for (const auto &entry : tree_to_center) { - const std::array ¢er = entry.second; + t8_forest_element_points_inside (forest, which_tree, element, data->points_flat, data->num_points, + points_inside.data (), data->tolerance); - const cell_coordinate cell = point_to_bbox_cell (center, bounds, level); + int count = 0; + t8_locidx_t last_inside = -1; - const auto key = cell_to_key (cell); - - const auto insert_result = occupied_cells.insert (key); - - if (!insert_result.second) { - return false; + for (int ipoint = 0; ipoint < data->num_points; ++ipoint) { + if (points_inside[static_cast (ipoint)]) { + count++; + last_inside = static_cast (ipoint); } } - return true; -} + int sum_of_elems = std::accumulate (points_inside.begin (), points_inside.end (), 0); -int -find_required_refinement_level (const std::map> &tree_to_center, - const double bounds[6], const int max_level) -{ - for (int level = 0; level <= max_level; ++level) { - if (refinement_level_is_unique (tree_to_center, bounds, level)) { - return level; - } + if (contained_tree_id != nullptr && count == 1) { + *contained_tree_id = last_inside; } - - return -1; + t8_productionf ("Points inside = %d\n", sum_of_elems); + return count; } -static t8_linearidx_t -cell_to_linear_id (const cell_coordinate &cell, const int level) +static int +count_tree_centers_inside_leaf2 (t8_forest_t forest, t8_locidx_t which_tree, const t8_element_t *element, + int num_points, double *points_flat) { - const t8_linearidx_t cells_per_direction = static_cast (1) << level; + std::vector points_inside (static_cast (num_points), 0); - return static_cast (cell.x) + cells_per_direction * static_cast (cell.y) - + cells_per_direction * cells_per_direction * static_cast (cell.z); + t8_forest_element_points_inside (forest, which_tree, element, points_flat, num_points, points_inside.data (), 1e-6); + + int sum_of_elems = std::accumulate (points_inside.begin (), points_inside.end (), 0); + t8_productionf ("Points inside = %d\n", sum_of_elems); + return sum_of_elems; +} +/* + * Adapt callback: + * + * Refine a bbox leaf if it contains more than one original tree center. + */ +static int +t8_adapt_refine ([[maybe_unused]] t8_forest_t forest, [[maybe_unused]] t8_forest_t forest_from, + [[maybe_unused]] const t8_locidx_t which_tree, [[maybe_unused]] const t8_eclass_t tree_class, + [[maybe_unused]] 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[]) +{ + // bbox_adapt_data *data = static_cast (t8_forest_get_user_data (forest_from)); + // if (data == nullptr) { + // data = static_cast (t8_forest_get_user_data (forest)); + // } + // T8_ASSERT (data != nullptr); + // + // const t8_element_t *element = elements[0]; + // const int level = scheme->element_get_level (tree_class, element); + // + // const int num_centers_inside = count_tree_centers_inside_leaf (forest_from, which_tree, element, data, nullptr); + // + // t8_productionf ("adapt check: bbox tree=%u, leaf id=%u, level=%i, centers_inside=%i%s\n", + // static_cast (which_tree), static_cast (lelement_id), level, num_centers_inside, + // num_centers_inside > 1 && level < data->max_level ? " -> refine" : " -> keep"); + // + // if (num_centers_inside > 1 && level < data->max_level) { + // data->refined_any = 1; + // return 1; // refine this bbox leaf + // } + // + // if (num_centers_inside > 1) { + // t8_productionf ("Warning: leaf contains %i centers, but max_level=%i was reached.\n", num_centers_inside, + // data->max_level); + // } + // + return 1; // keep this bbox leaf } std::map @@ -203,155 +227,246 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) { std::map tree_reindex; + t8_productionf ("starting tree reindexing\n"); + /* * 1. Get local bounding box of the original cmesh. */ double bounding_box[6]; t8_cmesh_get_local_bounding_box (cmesh, bounding_box); + //T8_ASSERT (t8_cmesh_get_local_bounding_box (cmesh, bounding_box)); + t8_productionf ("received local cmesh bounding box\n"); + t8_productionf ("physical bbox bounds: x=[%.17g, %.17g], y=[%.17g, %.17g], z=[%.17g, %.17g]\n", bounding_box[0], + bounding_box[1], bounding_box[2], bounding_box[3], bounding_box[4], bounding_box[5]); + + //const int use_quad_bbox = bbox_is_effectively_2d (bounding_box); /* * 2. Build auxiliary bounding-box cmesh. + * + * Important: The helper cmesh is a unit quad/hex. Physical tree centers from + * the original cmesh are normalized into this unit bbox before calling + * t8_forest_element_points_inside. This keeps the t8code inside-element + * implementation and avoids mixing physical coordinates from the original + * geometry with reference coordinates of the helper forest. */ - const auto vertices = bbox_bounds_to_hex_vertices (bounding_box); - t8_cmesh_t bbox_cmesh; t8_cmesh_init (&bbox_cmesh); + //if (use_quad_bbox) { + // const auto vertices = unit_quad_vertices (); + // t8_cmesh_set_tree_class (bbox_cmesh, 0, T8_ECLASS_QUAD); + // t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 4); + // t8_productionf ("initialized auxiliary bbox cmesh as unit QUAD\n"); + //} + //else { + // const auto vertices = unit_hex_vertices (); + // t8_cmesh_set_tree_class (bbox_cmesh, 0, T8_ECLASS_HEX); + // t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 8); + // t8_productionf ("initialized auxiliary bbox cmesh as unit HEX\n"); + //} + + std::vector vertices + = { bounding_box[0], bounding_box[2], bounding_box[4], bounding_box[1], bounding_box[3], bounding_box[5] }; + t8_cmesh_set_tree_class (bbox_cmesh, 0, T8_ECLASS_HEX); - t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 8); + t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 4); + t8_cmesh_register_geometry (bbox_cmesh); - sc_MPI_Comm bbox_comm = comm; + /* The helper bbox forest is local. Do not distribute this one helper tree over comm. */ + sc_MPI_Comm bbox_comm = sc_MPI_COMM_SELF; t8_cmesh_commit (bbox_cmesh, bbox_comm); + t8_productionf ("Committed auxiliary bbox cmesh\n"); + /* - * 3. Build level-0 forest from original cmesh to evaluate tree centers. + * 3. Build a level-0 forest from the original cmesh to evaluate tree centers. * * t8_forest_new_uniform takes ownership of a cmesh reference. - * Since cmesh belongs to the caller, we increase the refcount first. + * Since cmesh belongs to the caller, increase the refcount first. */ - t8_cmesh_ref (cmesh); + //t8_cmesh_ref (cmesh); + + t8_forest_t original_cmesh_forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default (), 0, 0, comm); - t8_forest_t center_forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default (), 0, 0, comm); + const t8_locidx_t num_cmesh_trees = t8_forest_get_num_global_trees (original_cmesh_forest); - const t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (center_forest); + t8_productionf ("created original level-0 forest with %u local tree(s)\n", static_cast (num_cmesh_trees)); /* * 4. Compute one physical center point per original local tree. */ std::map> tree_to_center; - for (t8_locidx_t itree = 0; itree < num_local_trees; ++itree) { - const std::array midpoint = compute_tree_midpoint_with_geometry (center_forest, itree); + for (t8_locidx_t itree = 0; itree < num_cmesh_trees; ++itree) { + const std::array midpoint = compute_tree_midpoint_with_geometry (original_cmesh_forest, itree); tree_to_center.insert ({ itree, midpoint }); - t8_productionf ("Tree Index: %u \t Tree Center: %f, %f, %f\n", static_cast (itree), midpoint[0], - midpoint[1], midpoint[2]); + t8_productionf ("Tree Index: %u \t physical Tree Center: %.17g, %.17g, %.17g\n", static_cast (itree), + midpoint[0], midpoint[1], midpoint[2]); } /* - * 5. Find the first bbox refinement level where each cell contains - * at most one tree center. + * 5. Normalize tree centers into the helper bbox coordinate system. + * + * Format: + * x0, y0, z0, x1, y1, z1, ... + * + * The index of a point in this array is also the original local tree id, + * because we fill it in tree-id order. */ - const int max_level = 29; + std::vector points_flat (static_cast (3 * num_cmesh_trees), 0.0); - const int required_level = find_required_refinement_level (tree_to_center, bounding_box, max_level); + const double tolerance = 1e-12; - if (required_level < 0) { - t8_productionf ("Could not find a refinement level where all tree centers are separated.\n" - "This can happen if two trees have exactly the same center.\n"); + for (const auto &entry : tree_to_center) { + const t8_locidx_t tree_id = entry.first; + const std::array &point = entry.second; - for (t8_locidx_t itree = 0; itree < num_local_trees; ++itree) { - tree_reindex[itree] = itree; - } + points_flat[static_cast (3 * tree_id + 0)] = point[0]; + points_flat[static_cast (3 * tree_id + 1)] = point[1]; + points_flat[static_cast (3 * tree_id + 2)] = point[2]; - t8_forest_unref (¢er_forest); + t8_productionf ("points_flat[%u] = %.17g, %.17g, %.17g\n", static_cast (tree_id), point[0], point[1], + point[2]); + } - /* - * bbox_cmesh was committed but not consumed by a forest in this branch. - */ - t8_cmesh_unref (&bbox_cmesh); + t8_productionf ("flattened %u tree center point(s)\n", static_cast (num_cmesh_trees)); - return tree_reindex; - } + /* + * 6. Create initial bbox forest at level 0. + * + * t8_forest_new_uniform takes ownership of bbox_cmesh and the scheme. + */ + t8_forest_t bbox_forest = t8_forest_new_uniform (bbox_cmesh, t8_scheme_new_default (), 0, 0, bbox_comm); - t8_productionf ("Required bounding-box refinement level: %i\n", required_level); + t8_productionf ("created initial level-0 bbox forest\n"); /* - * 6. Map each original tree center to the logical bbox-cell linear id - * at required_level. + * 7. Adaptively refine the bbox forest. * - * Since required_level passed the uniqueness test, each logical cell - * should contain at most one tree center. + * A leaf is refined iff it contains more than one normalized tree center. */ - std::map linear_id_to_tree; + bbox_adapt_data adapt_data; + adapt_data.points_flat = points_flat.data (); + adapt_data.num_points = static_cast (num_cmesh_trees); + adapt_data.tolerance = tolerance; + adapt_data.max_level = 29; + //adapt_data.refined_any = 0; + + int refinement_pass = 0; + + int num_bbox_trees = t8_forest_get_num_global_trees (bbox_forest); + + t8_productionf ("Now using own implementation of pointsinside\n"); + for (t8_locidx_t itree = 0; itree < num_bbox_trees; itree++) { + int num_elements_per_tree = t8_forest_get_tree_num_leaf_elements (bbox_forest, itree); + for (int ielement = 0; ielement < num_elements_per_tree; ielement++) { + const t8_element_t *element = t8_forest_get_leaf_element_in_tree (bbox_forest, itree, ielement); + count_tree_centers_inside_leaf2 (bbox_forest, itree, element, num_cmesh_trees, points_flat.data ()); + } + } - for (const auto &entry : tree_to_center) { - const t8_locidx_t original_tree = entry.first; - const std::array ¢er = entry.second; + while (true) { - const cell_coordinate cell = point_to_bbox_cell (center, bounding_box, required_level); + const t8_locidx_t leaf_count_before = t8_forest_get_local_num_leaf_elements (bbox_forest); - const t8_linearidx_t linear_id = cell_to_linear_id (cell, required_level); + t8_productionf ("starting refinement pass %i\n", refinement_pass); - const auto insert_result = linear_id_to_tree.insert ({ linear_id, original_tree }); + /* The callback reads user data from forest_from, matching normal t8code adapt usage. */ + t8_forest_set_user_data (bbox_forest, &adapt_data); - if (!insert_result.second) { - /* - * This should not happen because required_level was checked. - */ - t8_productionf ("Unexpected collision: more than one center in cell (%u, %u, %u), linear id %li.\n", - static_cast (cell.x), static_cast (cell.y), static_cast (cell.z), - static_cast (linear_id)); - } + t8_forest_t adapted_bbox_forest = t8_forest_new_adapt (bbox_forest, t8_adapt_refine, + 0, // non-recursive: one refinement step per loop iteration + 0, // no face ghosts + &adapt_data); - t8_productionf ("Tree %u lies in bbox cell (%u, %u, %u), linear id %li at level %i\n", - static_cast (original_tree), static_cast (cell.x), - static_cast (cell.y), static_cast (cell.z), static_cast (linear_id), - required_level); + const t8_locidx_t leaf_count_after = t8_forest_get_local_num_leaf_elements (adapted_bbox_forest); + + t8_productionf ("refinement pass %i leaf count: before=%u, after=%u\n", refinement_pass, + static_cast (leaf_count_before), static_cast (leaf_count_after)); + + /* + * We are done if no leaf requested refinement in this pass. + */ + //if (!adapt_data.refined_any) { + // t8_productionf ("refinement pass %i finished without further refinement; stopping\n", refinement_pass); + // + // t8_forest_unref (&bbox_forest); + // bbox_forest = adapted_bbox_forest; + // break; + //} + + t8_productionf ("refinement pass %i refined at least one leaf; continuing\n", refinement_pass); + + t8_forest_unref (&bbox_forest); + bbox_forest = adapted_bbox_forest; + + ++refinement_pass; } - /* - * 7. Create the refined bbox forest. - * - * This forest is only used to get the SFC order of the bbox leaves. - * t8_forest_new_uniform takes ownership of bbox_cmesh and the scheme. - */ - t8_forest_t bbox_forest = t8_forest_new_uniform (bbox_cmesh, t8_scheme_new_default (), required_level, 0, bbox_comm); + t8_productionf ("adaptive refinement finished after %i pass(es)\n", refinement_pass + 1); /* - * 8. Iterate over bbox leaves in t8code order. + * 8. Iterate over final bbox leaves in t8code/SFC order. * - * We do not evaluate geometry on bbox_forest here. - * Instead, we use each leaf's linear id at required_level. + * Since the adaptive refinement loop stopped, every leaf should contain + * at most one original tree center. Occupied leaves receive the next SFC + * index. */ + bbox_adapt_data final_check_data; + final_check_data.points_flat = points_flat.data (); + final_check_data.num_points = static_cast (num_cmesh_trees); + final_check_data.tolerance = tolerance; + final_check_data.max_level = 29; + //final_check_data.refined_any = 0; + t8_locidx_t next_sfc_index = 0; const t8_locidx_t num_bbox_local_trees = t8_forest_get_num_local_trees (bbox_forest); + t8_productionf ("iterating final bbox forest with %u local bbox tree(s)\n", + static_cast (num_bbox_local_trees)); + for (t8_locidx_t bbox_tree = 0; bbox_tree < num_bbox_local_trees; ++bbox_tree) { t8_element_array_t *leaf_elements = t8_forest_tree_get_leaf_elements (bbox_forest, bbox_tree); + t8_productionf ("iterating leaves of bbox tree %u\n", static_cast (bbox_tree)); for (auto leaf_it = t8_element_array_begin (leaf_elements); leaf_it != t8_element_array_end (leaf_elements); ++leaf_it) { - const t8_linearidx_t leaf_linear_id = leaf_it.get_linear_id_at_level (required_level); + const t8_element_t *leaf = *leaf_it; - const auto found_tree = linear_id_to_tree.find (leaf_linear_id); + t8_locidx_t contained_tree = -1; - if (found_tree == linear_id_to_tree.end ()) { + const int num_centers_inside + = count_tree_centers_inside_leaf1 (bbox_forest, bbox_tree, leaf, &final_check_data, &contained_tree); + + if (num_centers_inside == 0) { + t8_productionf ("final leaf in bbox tree %u is empty; skipping\n", static_cast (bbox_tree)); + continue; + } + + t8_productionf ("final leaf in bbox tree %u contains %i center(s); contained_tree=%i\n", + static_cast (bbox_tree), num_centers_inside, static_cast (contained_tree)); + + if (num_centers_inside > 1) { /* - * Empty bbox leaf. No original tree center lies here. + * This should only happen if max_level was reached, if two centers are + * identical, or if centers lie exactly on refinement boundaries and are + * considered inside multiple elements by the t8code predicate. */ + t8_productionf ("Warning: final bbox leaf still contains %i tree centers. " + "Using first detected tree only would be unsafe.\n", + num_centers_inside); + continue; } - const t8_locidx_t original_tree = found_tree->second; - - tree_reindex[original_tree] = next_sfc_index; + tree_reindex[contained_tree] = next_sfc_index; - t8_productionf ("bbox leaf linear id %li maps original tree %u to new SFC index %u\n", - static_cast (leaf_linear_id), static_cast (original_tree), + t8_productionf ("Original tree %u -> new SFC index %u\n", static_cast (contained_tree), static_cast (next_sfc_index)); ++next_sfc_index; @@ -361,24 +476,22 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) /* * 9. Sanity check. */ - if (tree_reindex.size () != static_cast (num_local_trees)) { + if (tree_reindex.size () != static_cast (num_cmesh_trees)) { t8_productionf ("Warning: only mapped %u of %u local trees.\n", static_cast (tree_reindex.size ()), - static_cast (num_local_trees)); + static_cast (num_cmesh_trees)); + } + else { + t8_productionf ("successfully mapped all %u local tree(s)\n", static_cast (num_cmesh_trees)); } /* * 10. Cleanup. */ + t8_productionf ("cleaning up bbox and original forests\n"); t8_forest_unref (&bbox_forest); - t8_forest_unref (¢er_forest); + t8_forest_unref (&original_cmesh_forest); + + t8_productionf ("tree reindexing finished\n"); return tree_reindex; } - -//void -//t8_perform_reindex (t8_cmesh_t cmesh, std::map tree_reindex_map) -//{ -// int num_trees = cmesh->num_trees; -// for (t8_locidx_t itree = 0; itree < num_trees; itree++) { -// } -//} From 54bfb091f3a7951e460a505a2cd50d6d383a657a Mon Sep 17 00:00:00 2001 From: Tu Nguyen Xuan Date: Mon, 6 Jul 2026 13:08:10 +0200 Subject: [PATCH 07/10] Reindex performed on stash --- .../t8_cmesh_tree_reindex.cxx | 484 +++++++++--------- test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx | 34 +- 2 files changed, 259 insertions(+), 259 deletions(-) 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 index 876db975ef..1f3fec8747 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -23,79 +23,25 @@ #include #include -#include #include +#include +#include +#include +#include #include +#include #include -#include -#include #include -#include -#include +#include +#include -#include #include -#include #include -#include #include - -//static int -//bbox_is_effectively_2d (const double bounds[6]) -//{ -// const double dz = std::fabs (bounds[5] - bounds[4]); -// const double scale = std::max ({ 1.0, std::fabs (bounds[4]), std::fabs (bounds[5]) }); -// return dz <= 1e-14 * scale; -//} - -//static double -//normalize_coordinate (const double value, const double min_value, const double max_value, const double tolerance) -//{ -// const double length = max_value - min_value; -// T8_ASSERT (std::fabs (length) > tolerance); -// -// double normalized = (value - min_value) / length; -// -// /* Avoid losing points due to tiny roundoff outside the unit bbox. */ -// if (normalized < 0.0 && normalized >= -tolerance) { -// normalized = 0.0; -// } -// if (normalized > 1.0 && normalized <= 1.0 + tolerance) { -// normalized = 1.0; -// } -// -// return normalized; -//} - -//static std::array -//unit_quad_vertices () -//{ -// return { -// 0.0, 0.0, 0.0, // vertex 0 -// 1.0, 0.0, 0.0, // vertex 1 -// 0.0, 1.0, 0.0, // vertex 2 -// 1.0, 1.0, 0.0 // vertex 3 -// }; -//} -// -//static std::array -//unit_hex_vertices () -//{ -// return { -// 0.0, 0.0, 0.0, // vertex 0 -// 1.0, 0.0, 0.0, // vertex 1 -// 0.0, 1.0, 0.0, // vertex 2 -// 1.0, 1.0, 0.0, // vertex 3 -// -// 0.0, 0.0, 1.0, // vertex 4 -// 1.0, 0.0, 1.0, // vertex 5 -// 0.0, 1.0, 1.0, // vertex 6 -// 1.0, 1.0, 1.0 // vertex 7 -// }; -//} +#include std::array -compute_tree_midpoint_with_geometry (t8_forest_t forest, t8_locidx_t local_tree_id) +compute_tree_midpoint_with_geometry1 (t8_forest_t forest, t8_locidx_t local_tree_id) { std::array physical_midpoint = { 0.0, 0.0, 0.0 }; @@ -109,123 +55,112 @@ compute_tree_midpoint_with_geometry (t8_forest_t forest, t8_locidx_t local_tree_ return physical_midpoint; } +std::array +compute_tree_midpoint_with_geometry (t8_forest_t forest, t8_locidx_t local_tree_id) +{ +} /* * Data passed to the adapt callback. * - * points_flat stores all normalized original tree centers as: + * points_flat stores all original tree centers in physical coordinates: * * x0, y0, z0, x1, y1, z1, ... * - * The helper bbox cmesh is deliberately built as a unit quad/hex. Therefore - * these points are normalized into [0, 1]^d before they are passed to - * t8_forest_element_points_inside. + * The point index is equal to the original local tree id. */ struct bbox_adapt_data { - const double *points_flat; + double *points_flat; int num_points; double tolerance; - int max_level; - /* Set to 1 by the callback whenever at least one leaf is refined. */ - //int refined_any; + /* + * Set to 1 by the callback whenever at least one leaf is refined. + */ + int refined_any; }; /* * Count how many original tree centers are inside a given bbox leaf. * - * This intentionally uses the t8code implementation t8_forest_element_points_inside. - * The points must be in the coordinate system of the helper bbox forest. Since - * the helper bbox forest is the unit quad/hex, points_flat contains normalized - * coordinates, not physical coordinates. - * * If contained_tree_id is not nullptr and exactly one point is inside, - * contained_tree_id will be set to the corresponding original tree id. + * contained_tree_id is set to the corresponding original local tree id. */ static int -count_tree_centers_inside_leaf1 (t8_forest_t forest, t8_locidx_t which_tree, const t8_element_t *element, - const bbox_adapt_data *data, t8_locidx_t *contained_tree_id) +count_tree_centers_inside_leaf (t8_forest_t forest, t8_locidx_t which_tree, const t8_element_t *element, int num_points, + double *points_flat, double tolerance, t8_locidx_t *contained_tree_id) { - std::vector points_inside (static_cast (data->num_points), 0); + std::vector points_inside (static_cast (num_points), 0); - t8_forest_element_points_inside (forest, which_tree, element, data->points_flat, data->num_points, - points_inside.data (), data->tolerance); + t8_forest_element_points_inside (forest, which_tree, element, points_flat, num_points, points_inside.data (), + tolerance); int count = 0; - t8_locidx_t last_inside = -1; + t8_locidx_t last_inside_tree = -1; - for (int ipoint = 0; ipoint < data->num_points; ++ipoint) { + for (int ipoint = 0; ipoint < num_points; ++ipoint) { if (points_inside[static_cast (ipoint)]) { - count++; - last_inside = static_cast (ipoint); + ++count; + last_inside_tree = static_cast (ipoint); } } - int sum_of_elems = std::accumulate (points_inside.begin (), points_inside.end (), 0); - - if (contained_tree_id != nullptr && count == 1) { - *contained_tree_id = last_inside; + if (count == 1 && contained_tree_id != nullptr) { + *contained_tree_id = last_inside_tree; } - t8_productionf ("Points inside = %d\n", sum_of_elems); + return count; } -static int -count_tree_centers_inside_leaf2 (t8_forest_t forest, t8_locidx_t which_tree, const t8_element_t *element, - int num_points, double *points_flat) -{ - std::vector points_inside (static_cast (num_points), 0); - - t8_forest_element_points_inside (forest, which_tree, element, points_flat, num_points, points_inside.data (), 1e-6); - - int sum_of_elems = std::accumulate (points_inside.begin (), points_inside.end (), 0); - t8_productionf ("Points inside = %d\n", sum_of_elems); - return sum_of_elems; -} /* - * Adapt callback: + * Adapt callback. * * Refine a bbox leaf if it contains more than one original tree center. */ static int -t8_adapt_refine ([[maybe_unused]] t8_forest_t forest, [[maybe_unused]] t8_forest_t forest_from, - [[maybe_unused]] const t8_locidx_t which_tree, [[maybe_unused]] const t8_eclass_t tree_class, - [[maybe_unused]] const t8_locidx_t lelement_id, [[maybe_unused]] const t8_scheme *scheme, +t8_adapt_refine ([[maybe_unused]] t8_forest_t forest, t8_forest_t forest_from, const t8_locidx_t which_tree, + const t8_eclass_t tree_class, [[maybe_unused]] 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[]) + t8_element_t *elements[]) { - // bbox_adapt_data *data = static_cast (t8_forest_get_user_data (forest_from)); - // if (data == nullptr) { - // data = static_cast (t8_forest_get_user_data (forest)); - // } - // T8_ASSERT (data != nullptr); - // - // const t8_element_t *element = elements[0]; - // const int level = scheme->element_get_level (tree_class, element); - // - // const int num_centers_inside = count_tree_centers_inside_leaf (forest_from, which_tree, element, data, nullptr); - // - // t8_productionf ("adapt check: bbox tree=%u, leaf id=%u, level=%i, centers_inside=%i%s\n", - // static_cast (which_tree), static_cast (lelement_id), level, num_centers_inside, - // num_centers_inside > 1 && level < data->max_level ? " -> refine" : " -> keep"); - // - // if (num_centers_inside > 1 && level < data->max_level) { - // data->refined_any = 1; - // return 1; // refine this bbox leaf - // } - // - // if (num_centers_inside > 1) { - // t8_productionf ("Warning: leaf contains %i centers, but max_level=%i was reached.\n", num_centers_inside, - // data->max_level); - // } - // - return 1; // keep this bbox leaf + bbox_adapt_data *data = static_cast (t8_forest_get_user_data (forest_from)); + + if (data == nullptr) { + data = static_cast (t8_forest_get_user_data (forest)); + } + + T8_ASSERT (data != nullptr); + + /* + * Important: + * Use the element passed by the callback. + * This element belongs to forest_from. + */ + const t8_element_t *element = elements[0]; + + const int level = scheme->element_get_level (tree_class, element); + + const int num_centers_inside = count_tree_centers_inside_leaf (forest_from, which_tree, element, data->num_points, + data->points_flat, data->tolerance, nullptr); + + const int should_refine = num_centers_inside > 1; + + t8_productionf ("adapt check: Tree=%u, Element=%u, level=%i, centers_inside=%i%s\n", + static_cast (which_tree), static_cast (lelement_id), level, num_centers_inside, + should_refine ? " -> refine" : " -> keep"); + + if (should_refine) { + data->refined_any = 1; + return 1; + } + + return 0; } -std::map +std::map t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) { - std::map tree_reindex; + std::map tree_reindex; t8_productionf ("starting tree reindexing\n"); @@ -234,68 +169,80 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) */ double bounding_box[6]; t8_cmesh_get_local_bounding_box (cmesh, bounding_box); - //T8_ASSERT (t8_cmesh_get_local_bounding_box (cmesh, bounding_box)); + t8_productionf ("received local cmesh bounding box\n"); t8_productionf ("physical bbox bounds: x=[%.17g, %.17g], y=[%.17g, %.17g], z=[%.17g, %.17g]\n", bounding_box[0], bounding_box[1], bounding_box[2], bounding_box[3], bounding_box[4], bounding_box[5]); - //const int use_quad_bbox = bbox_is_effectively_2d (bounding_box); - /* - * 2. Build auxiliary bounding-box cmesh. + * 2. Build auxiliary bbox cmesh. + * + * t8_geometry_linear_axis_aligned uses two vertices: * - * Important: The helper cmesh is a unit quad/hex. Physical tree centers from - * the original cmesh are normalized into this unit bbox before calling - * t8_forest_element_points_inside. This keeps the t8code inside-element - * implementation and avoids mixing physical coordinates from the original - * geometry with reference coordinates of the helper forest. + * min corner: x_min, y_min, z_min + * max corner: x_max, y_max, z_max */ t8_cmesh_t bbox_cmesh; t8_cmesh_init (&bbox_cmesh); - //if (use_quad_bbox) { - // const auto vertices = unit_quad_vertices (); - // t8_cmesh_set_tree_class (bbox_cmesh, 0, T8_ECLASS_QUAD); - // t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 4); - // t8_productionf ("initialized auxiliary bbox cmesh as unit QUAD\n"); - //} - //else { - // const auto vertices = unit_hex_vertices (); - // t8_cmesh_set_tree_class (bbox_cmesh, 0, T8_ECLASS_HEX); - // t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 8); - // t8_productionf ("initialized auxiliary bbox cmesh as unit HEX\n"); - //} - std::vector vertices = { bounding_box[0], bounding_box[2], bounding_box[4], bounding_box[1], bounding_box[3], bounding_box[5] }; - t8_cmesh_set_tree_class (bbox_cmesh, 0, T8_ECLASS_HEX); - t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 4); - t8_cmesh_register_geometry (bbox_cmesh); + const double dx = bounding_box[1] - bounding_box[0]; + const double dy = bounding_box[3] - bounding_box[2]; + const double dz = bounding_box[5] - bounding_box[4]; + + const double eps = 1e-14; + + const int active_dims = (std::abs (dx) > eps) + (std::abs (dy) > eps) + (std::abs (dz) > eps); + + t8_eclass_t bbox_eclass; + + switch (active_dims) { + case 3: + bbox_eclass = T8_ECLASS_HEX; + break; - /* The helper bbox forest is local. Do not distribute this one helper tree over comm. */ - sc_MPI_Comm bbox_comm = sc_MPI_COMM_SELF; + case 2: + bbox_eclass = T8_ECLASS_QUAD; + break; - t8_cmesh_commit (bbox_cmesh, bbox_comm); + case 1: + bbox_eclass = T8_ECLASS_LINE; + break; + + default: + SC_ABORT ("Bounding box has zero extent in all directions.\n"); + } + + t8_cmesh_set_tree_class (bbox_cmesh, 0, bbox_eclass); + t8_cmesh_set_tree_vertices (bbox_cmesh, 0, vertices.data (), 2); + t8_cmesh_register_geometry (bbox_cmesh); + + /* + * The bbox cmesh contains one helper tree. + */ + t8_cmesh_commit (bbox_cmesh, comm); t8_productionf ("Committed auxiliary bbox cmesh\n"); /* - * 3. Build a level-0 forest from the original cmesh to evaluate tree centers. + * 3. Build a level-0 forest from the original cmesh to compute + * one centroid per original tree. * * t8_forest_new_uniform takes ownership of a cmesh reference. * Since cmesh belongs to the caller, increase the refcount first. */ - //t8_cmesh_ref (cmesh); + t8_cmesh_ref (cmesh); t8_forest_t original_cmesh_forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default (), 0, 0, comm); - const t8_locidx_t num_cmesh_trees = t8_forest_get_num_global_trees (original_cmesh_forest); + const t8_locidx_t num_cmesh_trees = t8_forest_get_num_local_trees (original_cmesh_forest); t8_productionf ("created original level-0 forest with %u local tree(s)\n", static_cast (num_cmesh_trees)); /* - * 4. Compute one physical center point per original local tree. + * 4. Store all tree centers. */ std::map> tree_to_center; @@ -309,18 +256,17 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) } /* - * 5. Normalize tree centers into the helper bbox coordinate system. + * 5. Flatten tree centers for t8_forest_element_points_inside. * * Format: + * * x0, y0, z0, x1, y1, z1, ... * - * The index of a point in this array is also the original local tree id, - * because we fill it in tree-id order. + * Since we fill the array in local tree-id order, the point index is also + * the original local tree id. */ std::vector points_flat (static_cast (3 * num_cmesh_trees), 0.0); - const double tolerance = 1e-12; - for (const auto &entry : tree_to_center) { const t8_locidx_t tree_id = entry.first; const std::array &point = entry.second; @@ -340,69 +286,51 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) * * t8_forest_new_uniform takes ownership of bbox_cmesh and the scheme. */ - t8_forest_t bbox_forest = t8_forest_new_uniform (bbox_cmesh, t8_scheme_new_default (), 0, 0, bbox_comm); + t8_forest_t bbox_forest = t8_forest_new_uniform (bbox_cmesh, t8_scheme_new_default (), 0, 0, comm); t8_productionf ("created initial level-0 bbox forest\n"); /* * 7. Adaptively refine the bbox forest. * - * A leaf is refined iff it contains more than one normalized tree center. + * A bbox leaf is refined iff it contains more than one tree center. */ bbox_adapt_data adapt_data; adapt_data.points_flat = points_flat.data (); adapt_data.num_points = static_cast (num_cmesh_trees); - adapt_data.tolerance = tolerance; - adapt_data.max_level = 29; - //adapt_data.refined_any = 0; + adapt_data.tolerance = T8_PRECISION_SQRT_EPS; + adapt_data.refined_any = 0; int refinement_pass = 0; - int num_bbox_trees = t8_forest_get_num_global_trees (bbox_forest); - - t8_productionf ("Now using own implementation of pointsinside\n"); - for (t8_locidx_t itree = 0; itree < num_bbox_trees; itree++) { - int num_elements_per_tree = t8_forest_get_tree_num_leaf_elements (bbox_forest, itree); - for (int ielement = 0; ielement < num_elements_per_tree; ielement++) { - const t8_element_t *element = t8_forest_get_leaf_element_in_tree (bbox_forest, itree, ielement); - count_tree_centers_inside_leaf2 (bbox_forest, itree, element, num_cmesh_trees, points_flat.data ()); - } - } + t8_productionf ("starting refinement pass %i\n", refinement_pass); while (true) { + adapt_data.refined_any = 0; const t8_locidx_t leaf_count_before = t8_forest_get_local_num_leaf_elements (bbox_forest); - t8_productionf ("starting refinement pass %i\n", refinement_pass); - - /* The callback reads user data from forest_from, matching normal t8code adapt usage. */ - t8_forest_set_user_data (bbox_forest, &adapt_data); - - t8_forest_t adapted_bbox_forest = t8_forest_new_adapt (bbox_forest, t8_adapt_refine, - 0, // non-recursive: one refinement step per loop iteration - 0, // no face ghosts - &adapt_data); - - const t8_locidx_t leaf_count_after = t8_forest_get_local_num_leaf_elements (adapted_bbox_forest); - - t8_productionf ("refinement pass %i leaf count: before=%u, after=%u\n", refinement_pass, - static_cast (leaf_count_before), static_cast (leaf_count_after)); + t8_forest_t adapted_forest = t8_forest_new_adapt (bbox_forest, t8_adapt_refine, + 0, // non-recursive: one refinement step per loop iteration + 0, // no face ghosts + &adapt_data); /* - * We are done if no leaf requested refinement in this pass. + * Do not unref bbox_forest here. + * t8_forest_new_adapt takes ownership of the source forest. */ - //if (!adapt_data.refined_any) { - // t8_productionf ("refinement pass %i finished without further refinement; stopping\n", refinement_pass); - // - // t8_forest_unref (&bbox_forest); - // bbox_forest = adapted_bbox_forest; - // break; - //} + bbox_forest = adapted_forest; - t8_productionf ("refinement pass %i refined at least one leaf; continuing\n", refinement_pass); + const t8_locidx_t leaf_count_after = t8_forest_get_local_num_leaf_elements (bbox_forest); - t8_forest_unref (&bbox_forest); - bbox_forest = adapted_bbox_forest; + t8_productionf ("refinement pass %i leaf count: before=%u, after=%u, refined_any=%i\n", refinement_pass, + static_cast (leaf_count_before), static_cast (leaf_count_after), + adapt_data.refined_any); + + if (!adapt_data.refined_any) { + t8_productionf ("refinement pass %i finished without further refinement; stopping\n", refinement_pass); + break; + } ++refinement_pass; } @@ -410,66 +338,58 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) t8_productionf ("adaptive refinement finished after %i pass(es)\n", refinement_pass + 1); /* - * 8. Iterate over final bbox leaves in t8code/SFC order. + * 8. Iterate over the final adapted bbox forest in t8code/SFC order. * - * Since the adaptive refinement loop stopped, every leaf should contain - * at most one original tree center. Occupied leaves receive the next SFC - * index. + * Every final bbox leaf should contain at most one original tree center. + * If a leaf contains exactly one center, we assign the next SFC index to + * that original tree. */ - bbox_adapt_data final_check_data; - final_check_data.points_flat = points_flat.data (); - final_check_data.num_points = static_cast (num_cmesh_trees); - final_check_data.tolerance = tolerance; - final_check_data.max_level = 29; - //final_check_data.refined_any = 0; + t8_locidx_t new_tree_index = 0; - t8_locidx_t next_sfc_index = 0; + std::vector tree_was_mapped (static_cast (num_cmesh_trees), 0); const t8_locidx_t num_bbox_local_trees = t8_forest_get_num_local_trees (bbox_forest); - t8_productionf ("iterating final bbox forest with %u local bbox tree(s)\n", - static_cast (num_bbox_local_trees)); - - for (t8_locidx_t bbox_tree = 0; bbox_tree < num_bbox_local_trees; ++bbox_tree) { - t8_element_array_t *leaf_elements = t8_forest_tree_get_leaf_elements (bbox_forest, bbox_tree); - t8_productionf ("iterating leaves of bbox tree %u\n", static_cast (bbox_tree)); - - for (auto leaf_it = t8_element_array_begin (leaf_elements); leaf_it != t8_element_array_end (leaf_elements); - ++leaf_it) { - const t8_element_t *leaf = *leaf_it; - + for (t8_locidx_t bbox_itree = 0; bbox_itree < num_bbox_local_trees; ++bbox_itree) { + int 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_element_t *element = t8_forest_get_leaf_element_in_tree (bbox_forest, bbox_itree, ielement); t8_locidx_t contained_tree = -1; const int num_centers_inside - = count_tree_centers_inside_leaf1 (bbox_forest, bbox_tree, leaf, &final_check_data, &contained_tree); + = count_tree_centers_inside_leaf (bbox_forest, bbox_itree, element, static_cast (num_cmesh_trees), + adapt_data.points_flat, adapt_data.tolerance, &contained_tree); if (num_centers_inside == 0) { - t8_productionf ("final leaf in bbox tree %u is empty; skipping\n", static_cast (bbox_tree)); continue; } - t8_productionf ("final leaf in bbox tree %u contains %i center(s); contained_tree=%i\n", - static_cast (bbox_tree), num_centers_inside, static_cast (contained_tree)); - if (num_centers_inside > 1) { - /* - * This should only happen if max_level was reached, if two centers are - * identical, or if centers lie exactly on refinement boundaries and are - * considered inside multiple elements by the t8code predicate. - */ t8_productionf ("Warning: final bbox leaf still contains %i tree centers. " - "Using first detected tree only would be unsafe.\n", + "Skipping this leaf.\n", num_centers_inside); + continue; + } + + if (contained_tree < 0 || contained_tree >= num_cmesh_trees) { + t8_productionf ("Warning: invalid contained tree id %i. Skipping.\n", static_cast (contained_tree)); + continue; + } + if (tree_was_mapped[static_cast (contained_tree)]) { + t8_productionf ("Warning: original tree %u was already mapped. " + "Skipping duplicate occurrence.\n", + static_cast (contained_tree)); continue; } - tree_reindex[contained_tree] = next_sfc_index; + tree_reindex[contained_tree] = new_tree_index; + tree_was_mapped[static_cast (contained_tree)] = 1; t8_productionf ("Original tree %u -> new SFC index %u\n", static_cast (contained_tree), - static_cast (next_sfc_index)); + static_cast (new_tree_index)); - ++next_sfc_index; + ++new_tree_index; } } @@ -488,6 +408,7 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) * 10. Cleanup. */ t8_productionf ("cleaning up bbox and original forests\n"); + t8_forest_unref (&bbox_forest); t8_forest_unref (&original_cmesh_forest); @@ -495,3 +416,68 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) return tree_reindex; } + +static t8_gloidx_t +t8_cmesh_tree_reindex_lookup (const std::map &tree_reindex, const t8_gloidx_t old_id) +{ + const auto it = tree_reindex.find (old_id); + + T8_ASSERT (it != tree_reindex.end ()); + + return it->second; +} + +void +t8_cmesh_tree_perform_reindex_inplace (t8_stash_t stash, const std::map &tree_reindex) +{ + T8_ASSERT (stash != nullptr); + + /* Reindex tree class entries. */ + for (size_t iclass = 0; iclass < stash->classes.elem_count; ++iclass) { + t8_stash_class_struct_t *sclass = static_cast (sc_array_index (&stash->classes, iclass)); + + sclass->id = t8_cmesh_tree_reindex_lookup (tree_reindex, sclass->id); + } + + /* Reindex tree attributes. */ + for (size_t iattr = 0; iattr < stash->attributes.elem_count; ++iattr) { + t8_stash_attribute_struct_t *attr + = static_cast (sc_array_index (&stash->attributes, iattr)); + + attr->id = t8_cmesh_tree_reindex_lookup (tree_reindex, attr->id); + } + + /* Reindex face connections. */ + 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)); + + const t8_gloidx_t old_id1 = join->id1; + const t8_gloidx_t old_id2 = join->id2; + + const t8_gloidx_t new_id1 = t8_cmesh_tree_reindex_lookup (tree_reindex, old_id1); + + const t8_gloidx_t new_id2 = t8_cmesh_tree_reindex_lookup (tree_reindex, old_id2); + + const int old_face1 = join->face1; + const int old_face2 = join->face2; + + if (new_id1 <= new_id2) { + join->id1 = new_id1; + join->id2 = new_id2; + join->face1 = old_face1; + join->face2 = old_face2; + } + else { + join->id1 = new_id2; + join->id2 = new_id1; + join->face1 = old_face2; + join->face2 = old_face1; + } + } + + /* Restore the order expected by the stash/commit routines. */ + t8_stash_class_sort (stash); + t8_stash_joinface_sort (stash); + t8_stash_attribute_sort (stash); +} diff --git a/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx b/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx index d68b7d520f..873f19c0fa 100644 --- a/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx +++ b/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx @@ -4,7 +4,6 @@ #include -#include #include #include @@ -17,7 +16,7 @@ struct t8_test_cmesh_tree_reindex: public testing::Test const int do_partition = 0; const int periodic = 0; - cmesh = t8_cmesh_new_hypercube (T8_ECLASS_TRIANGLE, sc_MPI_COMM_SELF, do_partition, periodic, 0); + cmesh = t8_cmesh_new_hypercube (T8_ECLASS_TET, sc_MPI_COMM_SELF, do_partition, periodic, 0); ASSERT_NE (cmesh, nullptr); } @@ -34,19 +33,26 @@ struct t8_test_cmesh_tree_reindex: public testing::Test t8_cmesh_t cmesh = nullptr; }; -TEST_F (t8_test_cmesh_tree_reindex, sanity_check) +TEST_F (t8_test_cmesh_tree_reindex, sanity_check_tet_hypercube) { t8_productionf ("Test started\n"); const t8_locidx_t num_local_trees = t8_cmesh_get_num_local_trees (cmesh); - EXPECT_EQ (num_local_trees, 2); + /* + * A tetrahedral hypercube consists of 6 tetrahedron trees. + */ + EXPECT_EQ (num_local_trees, 6); std::map reindex = t8_cmesh_reindex_tree (cmesh, sc_MPI_COMM_SELF); + /* + * Every original local tree should receive exactly one new SFC index. + */ EXPECT_EQ (reindex.size (), static_cast (num_local_trees)); - std::set new_indices; + std::set old_tree_ids; + std::set new_sfc_indices; for (const auto &entry : reindex) { const t8_locidx_t old_tree_id = entry.first; @@ -58,17 +64,25 @@ TEST_F (t8_test_cmesh_tree_reindex, sanity_check) EXPECT_GE (new_sfc_index, 0); EXPECT_LT (new_sfc_index, num_local_trees); - new_indices.insert (new_sfc_index); + old_tree_ids.insert (old_tree_id); + new_sfc_indices.insert (new_sfc_index); t8_productionf ("old local tree id %li -> new local SFC index %li\n", static_cast (old_tree_id), static_cast (new_sfc_index)); } - EXPECT_EQ (new_indices.size (), static_cast (num_local_trees)); + /* + * The old tree ids should be exactly {0, 1, 2, 3, 4, 5}. + */ + EXPECT_EQ (old_tree_ids.size (), static_cast (num_local_trees)); /* - * Optional: for two trees, the SFC indices should be exactly {0, 1}. + * The new SFC indices should also be exactly {0, 1, 2, 3, 4, 5}. */ - EXPECT_TRUE (new_indices.count (0) == 1); - EXPECT_TRUE (new_indices.count (1) == 1); + EXPECT_EQ (new_sfc_indices.size (), static_cast (num_local_trees)); + + for (t8_locidx_t index = 0; index < num_local_trees; ++index) { + EXPECT_TRUE (old_tree_ids.count (index) == 1); + EXPECT_TRUE (new_sfc_indices.count (index) == 1); + } } From 03d22ae1ee8d4a9aa5d6b1bb37ccedd6f944ed2c Mon Sep 17 00:00:00 2001 From: Tu Nguyen Xuan Date: Thu, 9 Jul 2026 10:17:24 +0200 Subject: [PATCH 08/10] Working test --- .../t8_cmesh_tree_reindex.cxx | 253 ++++++++++++++---- .../t8_cmesh_tree_reindex.hxx | 2 +- test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx | 79 ++++-- 3 files changed, 265 insertions(+), 69 deletions(-) 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 index 1f3fec8747..11f1cf494b 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -34,30 +34,99 @@ #include #include #include +#include +#include #include #include #include #include - -std::array -compute_tree_midpoint_with_geometry1 (t8_forest_t forest, t8_locidx_t local_tree_id) +#include + +//std::array +//compute_tree_midpoint_with_geometry1 (t8_forest_t forest, t8_locidx_t local_tree_id) +//{ +// std::array physical_midpoint = { 0.0, 0.0, 0.0 }; +// +// const t8_element_t *element = t8_forest_get_leaf_element (forest, local_tree_id, 0); +// +// t8_forest_element_centroid (forest, local_tree_id, element, physical_midpoint.data ()); +// +// t8_productionf ("Computed centroid for local tree %u: %.17g, %.17g, %.17g\n", static_cast (local_tree_id), +// physical_midpoint[0], physical_midpoint[1], physical_midpoint[2]); +// +// return physical_midpoint; +//} + +static t8_eclass_t +get_tree_eclass_from_stash (const t8_stash_t stash, const t8_gloidx_t global_tree_id) { - std::array physical_midpoint = { 0.0, 0.0, 0.0 }; - - const t8_element_t *element = t8_forest_get_leaf_element (forest, local_tree_id, 0); + T8_ASSERT (stash != nullptr); - t8_forest_element_centroid (forest, local_tree_id, element, physical_midpoint.data ()); + for (size_t iclass = 0; iclass < stash->classes.elem_count; ++iclass) { + const t8_stash_class_struct_t *sclass + = static_cast (sc_array_index (&stash->classes, iclass)); - t8_productionf ("Computed centroid for local tree %u: %.17g, %.17g, %.17g\n", static_cast (local_tree_id), - physical_midpoint[0], physical_midpoint[1], physical_midpoint[2]); + if (sclass->id == global_tree_id) { + return sclass->eclass; + } + } - return physical_midpoint; + SC_ABORTF ("Could not find eclass for global tree %lli in stash.\n", static_cast (global_tree_id)); } -std::array -compute_tree_midpoint_with_geometry (t8_forest_t forest, t8_locidx_t local_tree_id) +static std::map> +compute_tree_centers_from_stash (const t8_stash_t stash) { + T8_ASSERT (stash != nullptr); + + std::map> tree_to_center; + + for (size_t iattr = 0; iattr < stash->attributes.elem_count; ++iattr) { + const t8_stash_attribute_struct_t *attr + = static_cast (sc_array_index (&stash->attributes, iattr)); + + if (attr->package_id != t8_get_package_id ()) { + continue; + } + + if (attr->key != T8_CMESH_VERTICES_ATTRIBUTE_KEY) { + continue; + } + + T8_ASSERT (attr->attr_data != nullptr); + T8_ASSERT (attr->attr_size % (3 * sizeof (double)) == 0); + + const t8_gloidx_t global_tree_id = attr->id; + + const t8_eclass_t eclass = get_tree_eclass_from_stash (stash, global_tree_id); + + const int expected_num_vertices = t8_eclass_num_vertices[eclass]; + + T8_ASSERT (expected_num_vertices > 0); + T8_ASSERT (stored_num_vertices >= expected_num_vertices); + + const double *vertices = static_cast (attr->attr_data); + + std::array center = { 0.0, 0.0, 0.0 }; + + for (int ivert = 0; ivert < expected_num_vertices; ++ivert) { + center[0] += vertices[3 * ivert + 0]; + center[1] += vertices[3 * ivert + 1]; + center[2] += vertices[3 * ivert + 2]; + } + + center[0] /= expected_num_vertices; + center[1] /= expected_num_vertices; + center[2] /= expected_num_vertices; + + tree_to_center.emplace (global_tree_id, center); + + t8_productionf ("Computed center for global tree %lli: %.17g, %.17g, %.17g\n", + static_cast (global_tree_id), center[0], center[1], center[2]); + } + + return tree_to_center; } /* * Data passed to the adapt callback. @@ -157,6 +226,85 @@ t8_adapt_refine ([[maybe_unused]] t8_forest_t forest, t8_forest_t forest_from, c return 0; } +// Helper Function to get an array with all vertices +struct t8_stash_vertex +{ + t8_gloidx_t tree_id; + int local_vertex_id; + std::array coords; +}; + +// Get flat vertices in the form x0, y0, z0, x1, y1, z1, ... xn, yn, zn of every tree +static std::vector +get_flat_vertices_from_stash (t8_stash_t stash) +{ + std::vector flat_vertices; + for (size_t iattr = 0; iattr < stash->attributes.elem_count; ++iattr) { + const t8_stash_attribute_struct_t *attr + = static_cast (sc_array_index (&stash->attributes, iattr)); + + if (attr->package_id != t8_get_package_id ()) { + continue; + } + + if (attr->key != T8_CMESH_VERTICES_ATTRIBUTE_KEY) { + continue; + } + + T8_ASSERT (attr->attr_data != nullptr); + T8_ASSERT (attr->attr_size % sizeof (double) == 0); + T8_ASSERT (attr->attr_size % (3 * sizeof (double)) == 0); + + const double *vertices = static_cast (attr->attr_data); + + const size_t num_doubles = attr->attr_size / sizeof (double); + + flat_vertices.insert (flat_vertices.end (), vertices, vertices + num_doubles); + } + return flat_vertices; +} + +/** + * Get min or max bound coordinates from a vertex list in the form of in the form x0, y0, z0, x1, y1, z1, ... xn, yn, zn, if which_bound == 0, then it returns the min bounds, so x_min, y_min, z_min + */ +static std::array +get_bounds_from_vertices (std::vector &vertices, int which_bound = 0) +{ + std::array return_coords; + if (!which_bound) { + std::array min_coords = { std::numeric_limits::max (), std::numeric_limits::max (), + std::numeric_limits::max () }; + + for (size_t i = 0; i < vertices.size (); i += 3) { + const double x = vertices[i + 0]; + const double y = vertices[i + 1]; + const double z = vertices[i + 2]; + + min_coords[0] = std::min (min_coords[0], x); + min_coords[1] = std::min (min_coords[1], y); + min_coords[2] = std::min (min_coords[2], z); + return_coords = min_coords; + } + } + else { + std::array max_coords = { std::numeric_limits::min (), std::numeric_limits::min (), + std::numeric_limits::min () }; + + for (size_t i = 0; i < vertices.size (); i += 3) { + const double x = vertices[i + 0]; + const double y = vertices[i + 1]; + const double z = vertices[i + 2]; + + max_coords[0] = std::max (max_coords[0], x); + max_coords[1] = std::max (max_coords[1], y); + max_coords[2] = std::max (max_coords[2], z); + return_coords = max_coords; + } + } + + return return_coords; +} + std::map t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) { @@ -164,6 +312,19 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) t8_productionf ("starting tree reindexing\n"); + /** + * Iterating through all vertices of the cmesh in each tree. Store global maximum and minimum coordinates to later create the bbox cmesh, and also store tree centers all in one iteration. + */ + + t8_stash_t original_cmesh_stash = cmesh->stash; + + //vector with all vertices, flattened out, remapping to its original trees will be done later + std::vector flat_vertices = get_flat_vertices_from_stash (original_cmesh_stash); + + //Min and Max Bounds for later bounding box cmesh + std::array max_coords = get_bounds_from_vertices (flat_vertices, 1); + std::array min_coords = get_bounds_from_vertices (flat_vertices); + /* * 1. Get local bounding box of the original cmesh. */ @@ -171,8 +332,8 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) t8_cmesh_get_local_bounding_box (cmesh, bounding_box); t8_productionf ("received local cmesh bounding box\n"); - t8_productionf ("physical bbox bounds: x=[%.17g, %.17g], y=[%.17g, %.17g], z=[%.17g, %.17g]\n", bounding_box[0], - bounding_box[1], bounding_box[2], bounding_box[3], bounding_box[4], bounding_box[5]); + t8_productionf ("physical bbox bounds: x=[%.17g, %.17g], y=[%.17g, %.17g], z=[%.17g, %.17g]\n", min_coords[0], + max_coords[0], min_coords[1], max_coords[1], min_coords[2], max_coords[2]); /* * 2. Build auxiliary bbox cmesh. @@ -186,16 +347,17 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) t8_cmesh_init (&bbox_cmesh); std::vector vertices - = { bounding_box[0], bounding_box[2], bounding_box[4], bounding_box[1], bounding_box[3], bounding_box[5] }; + = { min_coords[0], min_coords[1], min_coords[2], max_coords[0], max_coords[1], max_coords[2] }; - const double dx = bounding_box[1] - bounding_box[0]; - const double dy = bounding_box[3] - bounding_box[2]; - const double dz = bounding_box[5] - bounding_box[4]; + const double dx = max_coords[0] - min_coords[0]; + const double dy = max_coords[1] - min_coords[1]; + const double dz = max_coords[2] - min_coords[2]; - const double eps = 1e-14; + const double eps = T8_PRECISION_SQRT_EPS; const int active_dims = (std::abs (dx) > eps) + (std::abs (dy) > eps) + (std::abs (dz) > eps); + t8_productionf ("Active Dimension %u\n", active_dims); t8_eclass_t bbox_eclass; switch (active_dims) { @@ -223,6 +385,7 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) * The bbox cmesh contains one helper tree. */ t8_cmesh_commit (bbox_cmesh, comm); + t8_cmesh_vtk_write_file (bbox_cmesh, "debug_bbox_cmesh"); t8_productionf ("Committed auxiliary bbox cmesh\n"); @@ -233,27 +396,13 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) * t8_forest_new_uniform takes ownership of a cmesh reference. * Since cmesh belongs to the caller, increase the refcount first. */ - t8_cmesh_ref (cmesh); - - t8_forest_t original_cmesh_forest = t8_forest_new_uniform (cmesh, t8_scheme_new_default (), 0, 0, comm); - - const t8_locidx_t num_cmesh_trees = t8_forest_get_num_local_trees (original_cmesh_forest); - - t8_productionf ("created original level-0 forest with %u local tree(s)\n", static_cast (num_cmesh_trees)); /* * 4. Store all tree centers. */ - std::map> tree_to_center; - - for (t8_locidx_t itree = 0; itree < num_cmesh_trees; ++itree) { - const std::array midpoint = compute_tree_midpoint_with_geometry (original_cmesh_forest, itree); - - tree_to_center.insert ({ itree, midpoint }); + std::map> tree_to_center = compute_tree_centers_from_stash (original_cmesh_stash); - t8_productionf ("Tree Index: %u \t physical Tree Center: %.17g, %.17g, %.17g\n", static_cast (itree), - midpoint[0], midpoint[1], midpoint[2]); - } + const t8_locidx_t num_cmesh_trees = static_cast (tree_to_center.size ()); /* * 5. Flatten tree centers for t8_forest_element_points_inside. @@ -265,18 +414,24 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) * Since we fill the array in local tree-id order, the point index is also * the original local tree id. */ - std::vector points_flat (static_cast (3 * num_cmesh_trees), 0.0); + std::vector points_flat; + points_flat.reserve (3 * tree_to_center.size ()); + + std::vector point_index_to_global_tree_id; + point_index_to_global_tree_id.reserve (tree_to_center.size ()); for (const auto &entry : tree_to_center) { - const t8_locidx_t tree_id = entry.first; + const t8_gloidx_t global_tree_id = entry.first; const std::array &point = entry.second; - points_flat[static_cast (3 * tree_id + 0)] = point[0]; - points_flat[static_cast (3 * tree_id + 1)] = point[1]; - points_flat[static_cast (3 * tree_id + 2)] = point[2]; + point_index_to_global_tree_id.push_back (global_tree_id); - t8_productionf ("points_flat[%u] = %.17g, %.17g, %.17g\n", static_cast (tree_id), point[0], point[1], - point[2]); + points_flat.push_back (point[0]); + points_flat.push_back (point[1]); + points_flat.push_back (point[2]); + + t8_productionf ("points_flat for global tree %lli = %.17g, %.17g, %.17g\n", static_cast (global_tree_id), + point[0], point[1], point[2]); } t8_productionf ("flattened %u tree center point(s)\n", static_cast (num_cmesh_trees)); @@ -334,6 +489,13 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) ++refinement_pass; } + t8_forest_vtk_write_file (bbox_forest, "debug_bbox_forest", 1, /* write_treeid */ + 1, /* write_mpirank */ + 1, /* write_level */ + 1, /* write_element_id */ + 0, /* write_ghosts */ + 0, /* num_data */ + nullptr); /* data */ t8_productionf ("adaptive refinement finished after %i pass(es)\n", refinement_pass + 1); @@ -383,7 +545,9 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) continue; } - tree_reindex[contained_tree] = new_tree_index; + const t8_gloidx_t old_global_tree_id = point_index_to_global_tree_id[static_cast (contained_tree)]; + + tree_reindex[old_global_tree_id] = static_cast (new_tree_index); tree_was_mapped[static_cast (contained_tree)] = 1; t8_productionf ("Original tree %u -> new SFC index %u\n", static_cast (contained_tree), @@ -407,10 +571,9 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) /* * 10. Cleanup. */ - t8_productionf ("cleaning up bbox and original forests\n"); + t8_productionf ("cleaning up bbox\n"); t8_forest_unref (&bbox_forest); - t8_forest_unref (&original_cmesh_forest); t8_productionf ("tree reindexing finished\n"); 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 index 89d345e016..54f2c2985c 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.hxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.hxx @@ -25,5 +25,5 @@ #include #include -std::map +std::map t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm); diff --git a/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx b/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx index 873f19c0fa..1048d662a7 100644 --- a/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx +++ b/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx @@ -1,24 +1,50 @@ #include -#include +#include #include #include +#include #include #include struct t8_test_cmesh_tree_reindex: public testing::Test { protected: + static void + add_tet_tree (t8_cmesh_t cmesh, const t8_gloidx_t global_tree_id, const std::array &vertices) + { + t8_cmesh_set_tree_class (cmesh, global_tree_id, T8_ECLASS_TET); + t8_cmesh_set_tree_vertices (cmesh, global_tree_id, vertices.data (), 4); + } + void SetUp () override { - const int do_partition = 0; - const int periodic = 0; - - cmesh = t8_cmesh_new_hypercube (T8_ECLASS_TET, sc_MPI_COMM_SELF, do_partition, periodic, 0); + t8_cmesh_init (&cmesh); ASSERT_NE (cmesh, nullptr); + + /* + * Build an uncommitted tetrahedral cmesh manually. + * + * Each tree has 4 vertices, stored as: + * + * x0, y0, z0, x1, y1, z1, ... + * + * The centers are distinct and lie inside the unit cube. + */ + add_tet_tree (cmesh, 0, { 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0 }); + + add_tet_tree (cmesh, 1, { 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 }); + + add_tet_tree (cmesh, 2, { 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 }); + + add_tet_tree (cmesh, 3, { 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 }); + + add_tet_tree (cmesh, 4, { 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 }); + + add_tet_tree (cmesh, 5, { 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0 }); } void @@ -33,55 +59,62 @@ struct t8_test_cmesh_tree_reindex: public testing::Test t8_cmesh_t cmesh = nullptr; }; -TEST_F (t8_test_cmesh_tree_reindex, sanity_check_tet_hypercube) +TEST_F (t8_test_cmesh_tree_reindex, sanity_check_uncommitted_tet_stash) { t8_productionf ("Test started\n"); - const t8_locidx_t num_local_trees = t8_cmesh_get_num_local_trees (cmesh); + constexpr t8_gloidx_t num_trees = 6; /* - * A tetrahedral hypercube consists of 6 tetrahedron trees. + * Important: + * The cmesh must still be uncommitted here. + * The reindexing function is supposed to read from cmesh->stash. */ - EXPECT_EQ (num_local_trees, 6); + EXPECT_FALSE (t8_cmesh_is_committed (cmesh)); - std::map reindex = t8_cmesh_reindex_tree (cmesh, sc_MPI_COMM_SELF); + std::map reindex = t8_cmesh_reindex_tree (cmesh, sc_MPI_COMM_SELF); + + /* + * Reindexing should not commit the original cmesh. + */ + EXPECT_FALSE (t8_cmesh_is_committed (cmesh)); /* - * Every original local tree should receive exactly one new SFC index. + * Every original global tree id should receive exactly one new SFC index. */ - EXPECT_EQ (reindex.size (), static_cast (num_local_trees)); + EXPECT_EQ (reindex.size (), static_cast (num_trees)); - std::set old_tree_ids; - std::set new_sfc_indices; + std::set old_tree_ids; + std::set new_sfc_indices; for (const auto &entry : reindex) { - const t8_locidx_t old_tree_id = entry.first; - const t8_locidx_t new_sfc_index = entry.second; + const t8_gloidx_t old_tree_id = entry.first; + const t8_gloidx_t new_sfc_index = entry.second; EXPECT_GE (old_tree_id, 0); - EXPECT_LT (old_tree_id, num_local_trees); + EXPECT_LT (old_tree_id, num_trees); EXPECT_GE (new_sfc_index, 0); - EXPECT_LT (new_sfc_index, num_local_trees); + EXPECT_LT (new_sfc_index, num_trees); old_tree_ids.insert (old_tree_id); new_sfc_indices.insert (new_sfc_index); - t8_productionf ("old local tree id %li -> new local SFC index %li\n", static_cast (old_tree_id), - static_cast (new_sfc_index)); + t8_productionf ("old global tree id %lli -> new global SFC index %lli\n", static_cast (old_tree_id), + static_cast (new_sfc_index)); } /* * The old tree ids should be exactly {0, 1, 2, 3, 4, 5}. */ - EXPECT_EQ (old_tree_ids.size (), static_cast (num_local_trees)); + EXPECT_EQ (old_tree_ids.size (), static_cast (num_trees)); /* * The new SFC indices should also be exactly {0, 1, 2, 3, 4, 5}. */ - EXPECT_EQ (new_sfc_indices.size (), static_cast (num_local_trees)); + EXPECT_EQ (new_sfc_indices.size (), static_cast (num_trees)); - for (t8_locidx_t index = 0; index < num_local_trees; ++index) { + for (t8_gloidx_t index = 0; index < num_trees; ++index) { EXPECT_TRUE (old_tree_ids.count (index) == 1); EXPECT_TRUE (new_sfc_indices.count (index) == 1); } From e68eb5331af8e9266281b4ea87bb7bef35a49354 Mon Sep 17 00:00:00 2001 From: Tu Nguyen Xuan Date: Thu, 9 Jul 2026 11:23:31 +0200 Subject: [PATCH 09/10] deleted vtk file writes --- src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx | 8 -------- 1 file changed, 8 deletions(-) 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 index 11f1cf494b..658517578e 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -385,7 +385,6 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) * The bbox cmesh contains one helper tree. */ t8_cmesh_commit (bbox_cmesh, comm); - t8_cmesh_vtk_write_file (bbox_cmesh, "debug_bbox_cmesh"); t8_productionf ("Committed auxiliary bbox cmesh\n"); @@ -489,13 +488,6 @@ t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) ++refinement_pass; } - t8_forest_vtk_write_file (bbox_forest, "debug_bbox_forest", 1, /* write_treeid */ - 1, /* write_mpirank */ - 1, /* write_level */ - 1, /* write_element_id */ - 0, /* write_ghosts */ - 0, /* num_data */ - nullptr); /* data */ t8_productionf ("adaptive refinement finished after %i pass(es)\n", refinement_pass + 1); From 25114ab08626e8a65cb7e058384bee664b7990a2 Mon Sep 17 00:00:00 2001 From: Tu Nguyen Xuan Date: Thu, 9 Jul 2026 11:24:57 +0200 Subject: [PATCH 10/10] removed old assertion --- src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx | 1 - 1 file changed, 1 deletion(-) 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 index 658517578e..5af3a6d5f5 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -104,7 +104,6 @@ compute_tree_centers_from_stash (const t8_stash_t stash) const int expected_num_vertices = t8_eclass_num_vertices[eclass]; T8_ASSERT (expected_num_vertices > 0); - T8_ASSERT (stored_num_vertices >= expected_num_vertices); const double *vertices = static_cast (attr->attr_data);