diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 6d5197d07e..7b61b3c0d5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -151,6 +151,7 @@ target_sources( T8 PRIVATE t8_cmesh/t8_cmesh_internal/t8_cmesh_partition.cxx t8_cmesh/t8_cmesh_internal/t8_cmesh_stash.c t8_cmesh/t8_cmesh_internal/t8_cmesh_trees.cxx + t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_conn_tree_to_vertex.cxx t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_conn_vertex_to_tree.cxx t8_cmesh/t8_cmesh_vertex_connectivity/t8_cmesh_vertex_connectivity.cxx diff --git a/src/t8_cmesh/t8_cmesh.cxx b/src/t8_cmesh/t8_cmesh.cxx index 27c7e7cd79..e0df232faa 100644 --- a/src/t8_cmesh/t8_cmesh.cxx +++ b/src/t8_cmesh/t8_cmesh.cxx @@ -152,6 +152,12 @@ t8_cmesh_disable_negative_volume_check ([[maybe_unused]] t8_cmesh_t cmesh) #endif } +void +t8_cmesh_enable_tree_reordering (t8_cmesh_t cmesh) +{ + cmesh->reindex_trees = 1; +} + #if T8_ENABLE_DEBUG int t8_cmesh_validate_geometry (const t8_cmesh_t cmesh, const int check_for_negative_volume) @@ -239,6 +245,7 @@ t8_cmesh_init (t8_cmesh_t *pcmesh) #if T8_ENABLE_DEBUG cmesh->negative_volume_check = 1; #endif /* T8_ENABLE_DEBUG */ + cmesh->reindex_trees = 0; T8_ASSERT (t8_cmesh_is_initialized (cmesh)); } @@ -2254,8 +2261,8 @@ t8_cmesh_uniform_bounds_from_partition (const t8_cmesh_t cmesh, const t8_gloidx_ #endif } /* End sending of end message */ - } /* End loop over processes */ - } /* if (num_pure_local_trees > 0) */ + } /* End loop over processes */ + } /* if (num_pure_local_trees > 0) */ if (this_proc_is_empty) { /* If this proc is empty, we only expect a start message. diff --git a/src/t8_cmesh/t8_cmesh.h b/src/t8_cmesh/t8_cmesh.h index 71c6520fe2..393711bc53 100644 --- a/src/t8_cmesh/t8_cmesh.h +++ b/src/t8_cmesh/t8_cmesh.h @@ -132,6 +132,13 @@ t8_cmesh_stash_is_empty (const t8_cmesh_t cmesh); void t8_cmesh_disable_negative_volume_check (t8_cmesh_t cmesh); +/** + * Enable localitly based indexing of trees during \ref t8_cmesh_commit. + * \param [in, out] cmesh + */ +void +t8_cmesh_enable_tree_reordering (t8_cmesh_t cmesh); + #if T8_ENABLE_DEBUG /** Check the geometry of the mesh for validity, this means checking if trees and their geometries * are compatible and if they have negative volume. @@ -140,7 +147,6 @@ t8_cmesh_disable_negative_volume_check (t8_cmesh_t cmesh); * \return True if the geometry of the cmesh is valid. */ int - t8_cmesh_validate_geometry (const t8_cmesh_t cmesh, const int check_for_negative_volume); #endif diff --git a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_commit.cxx b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_commit.cxx index 7d0f789302..5c48d7f346 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_commit.cxx +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_commit.cxx @@ -36,6 +36,8 @@ #include #include #include +#include +#include /** * A struct to hold the information about a ghost facejoin. @@ -148,6 +150,16 @@ t8_cmesh_commit_replicated_new (t8_cmesh_t cmesh) t8_stash_class_struct_t *entry; t8_locidx_t num_trees = class_entries->elem_count, ltree; + if (cmesh->reindex_trees) { + auto start = std::chrono::high_resolution_clock::now (); + std::map reindexing_map = t8_cmesh_reindex_tree (cmesh); + + t8_cmesh_tree_perform_reindex_inplace (stash, reindexing_map); + auto stop = std::chrono::high_resolution_clock::now (); + auto duration = duration_cast (stop - start); + t8_productionf ("Time taken: %lims\n", duration.count ()); + } + t8_cmesh_trees_init (&cmesh->trees, 1, num_trees, 0); t8_cmesh_trees_start_part (cmesh->trees, 0, 0, num_trees, 0, 0, 1); /* set tree classes */ diff --git a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx new file mode 100644 index 0000000000..b0f2148a49 --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.cxx @@ -0,0 +1,512 @@ +/* + 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +//static t8_eclass_t +//get_tree_eclass_from_stash (const t8_stash_t stash, const t8_gloidx_t global_tree_id) +//{ +// T8_ASSERT (stash != nullptr); +// +// 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)); +// +// if (sclass->id == global_tree_id) { +// return sclass->eclass; +// } +// } +// +// SC_ABORTF ("Could not find eclass for global tree %lli in stash.\n", static_cast (global_tree_id)); +//} + +// the data for each element. Each element has a list of particles it contains +struct element_data +{ + std::vector> midpoints; +}; + +// the forest data contains the data of the elements and a flag, which lets us know when we are finished +// in the transfer_points function you can check, if one of the new elements contains more than one particle +// and then you set finished to false. At the end of the loop you check the state of finished +struct forest_data +{ + element_data *elem_data; + t8_locidx_t num_elements; + int finished; +}; + +static forest_data * +forest_data_new (const t8_locidx_t num_elements) +{ + forest_data *data = T8_ALLOC (forest_data, 1); + + data->num_elements = num_elements; + data->finished = 1; + data->elem_data = T8_ALLOC (element_data, num_elements); + + for (t8_locidx_t ielement = 0; ielement < num_elements; ++ielement) { + new (&data->elem_data[ielement]) element_data (); + } + + return data; +} + +static void +forest_data_destroy (forest_data *data) +{ + if (data == nullptr) { + return; + } + + for (t8_locidx_t ielement = 0; ielement < data->num_elements; ++ielement) { + data->elem_data[ielement].~element_data (); + } + + T8_FREE (data->elem_data); + T8_FREE (data); +} + +static int +t8_adapt_refine ([[maybe_unused]] t8_forest_t forest, t8_forest_t forest_from, const t8_locidx_t which_tree, + [[maybe_unused]] const t8_eclass_t tree_class, [[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[]) +{ + forest_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); + t8_locidx_t element_index = t8_forest_get_tree_element_offset (forest_from, which_tree) + lelement_id; + if (data->elem_data[element_index].midpoints.size () > 1) { + return 1; + } + + return 0; +} + +static void +transfer_points (t8_forest_t forest_old, t8_forest_t forest_new, const t8_locidx_t which_tree, + [[maybe_unused]] const t8_eclass_t tree_class, [[maybe_unused]] const t8_scheme *scheme, + const int refine, [[maybe_unused]] const int num_old_elements, const t8_locidx_t first_old_element, + const int num_new_elements, const t8_locidx_t first_new_element) +{ + const forest_data *data_old = static_cast (t8_forest_get_user_data (forest_old)); + forest_data *data_new = static_cast (t8_forest_get_user_data (forest_new)); + + if (refine == 0) { + + t8_locidx_t old_index = t8_forest_get_tree_element_offset (forest_old, which_tree) + first_old_element; + t8_locidx_t new_index = t8_forest_get_tree_element_offset (forest_new, which_tree) + first_new_element; + + data_new->elem_data[new_index] = data_old->elem_data[old_index]; + + if (data_new->elem_data[new_index].midpoints.size () > 1) { + data_new->finished = 0; + } + } + + else { + + t8_locidx_t old_index = t8_forest_get_tree_element_offset (forest_old, which_tree) + first_old_element; + const element_data &old_element_data = data_old->elem_data[old_index]; + + if (old_element_data.midpoints.empty ()) { + return; + } + + std::vector point_coords (3 * old_element_data.midpoints.size ()); + std::vector point_was_copied (old_element_data.midpoints.size (), 0); + + for (std::size_t ipoint = 0; ipoint < old_element_data.midpoints.size (); ++ipoint) { + const t8_3D_vec &point = old_element_data.midpoints[ipoint].second; + + point_coords[3 * ipoint + 0] = point[0]; + point_coords[3 * ipoint + 1] = point[1]; + point_coords[3 * ipoint + 2] = point[2]; + } + + for (int inew = 0; inew < num_new_elements; ++inew) { + t8_locidx_t new_tree_leaf_index = first_new_element + inew; + t8_locidx_t new_index = t8_forest_get_tree_element_offset (forest_new, which_tree) + new_tree_leaf_index; + + const t8_element_t *new_element + = t8_forest_get_leaf_element_in_tree (forest_new, which_tree, new_tree_leaf_index); + + std::vector point_inside (old_element_data.midpoints.size (), 0); + + t8_forest_element_points_inside (forest_new, which_tree, new_element, point_coords.data (), + static_cast (old_element_data.midpoints.size ()), point_inside.data (), + 0.0f); + + for (std::size_t ipoint = 0; ipoint < old_element_data.midpoints.size (); ++ipoint) { + if (point_inside[ipoint] && !point_was_copied[ipoint]) { + data_new->elem_data[new_index].midpoints.push_back (old_element_data.midpoints[ipoint]); + point_was_copied[ipoint] = 1; + } + } + + if (data_new->elem_data[new_index].midpoints.size () > 1) { + data_new->finished = false; + } + } + } +} + +std::map +t8_cmesh_reindex_tree (t8_cmesh_t cmesh, sc_MPI_Comm comm) +{ + std::map tree_reindex; + std::vector flat_vertices; + + t8_productionf ("starting tree reindexing\n"); + t8_stash_t original_cmesh_stash = cmesh->stash; + /** + * 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_3D_vec min_coords + = { std::numeric_limits::max (), std::numeric_limits::max (), std::numeric_limits::max () }; + + t8_3D_vec max_coords = { std::numeric_limits::lowest (), std::numeric_limits::lowest (), + std::numeric_limits::lowest () }; + + std::map tree_to_center; + + for (size_t iattr = 0; iattr < original_cmesh_stash->attributes.elem_count; ++iattr) { + const t8_stash_attribute_struct_t *attr + = static_cast (sc_array_index (&original_cmesh_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 int expected_num_vertices = attr->attr_size / (3 * sizeof (double)); + const std::span vertices (static_cast (attr->attr_data), expected_num_vertices); + + t8_3D_vec center = { 0.0, 0.0, 0.0 }; + + for (const auto &ivert : vertices) { + + t8_productionf ("Coordinates are (%f, %f, %f)\n", ivert[0], ivert[1], ivert[2]); + for (int idim = 0; idim < 3; ++idim) { + min_coords[idim] = std::min (min_coords[idim], ivert[idim]); + max_coords[idim] = std::max (max_coords[idim], ivert[idim]); + center[idim] += ivert[idim]; + } + } + + center[0] /= expected_num_vertices; + center[1] /= expected_num_vertices; + center[2] /= expected_num_vertices; + + tree_to_center.emplace (attr->id, center); + + t8_productionf ("Computed center for global tree %li: %f, %f, %.f\n", attr->id, center[0], center[1], center[2]); + } + + t8_productionf ("received local cmesh bounding box\n"); + 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]); + + t8_cmesh_t bbox_cmesh; + t8_cmesh_init (&bbox_cmesh); + + std::vector vertices + = { min_coords[0], min_coords[1], min_coords[2], max_coords[0], max_coords[1], max_coords[2] }; + + 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 = 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) { + case 3: + bbox_eclass = T8_ECLASS_HEX; + break; + + case 2: + bbox_eclass = T8_ECLASS_QUAD; + break; + + case 1: + bbox_eclass = T8_ECLASS_LINE; + break; + + default: + SC_ABORT ("Bounding box has zero extent in all directions.\n"); + } + + 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); + + t8_cmesh_commit (bbox_cmesh, comm); + + t8_productionf ("Committed auxiliary bbox cmesh\n"); + + const int num_cmesh_trees = tree_to_center.size (); + + t8_productionf ("flattened %u tree center point(s)\n", static_cast (num_cmesh_trees)); + + 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"); + + if (!t8_forest_write_vtk (bbox_forest, "bounding_box")) { + t8_productionf ("Could not write VTK file for forest"); + }; + + forest_data *data = forest_data_new (t8_forest_get_local_num_leaf_elements (bbox_forest)); + + T8_ASSERT (data->num_elements == 1); + + for (const auto &entry : tree_to_center) { + const t8_gloidx_t global_tree_id = entry.first; + const t8_3D_vec &midpoint = entry.second; + + data->elem_data[0].midpoints.emplace_back (std::make_pair (global_tree_id, midpoint)); + } + + if (data->elem_data[0].midpoints.size () > 1) { + data->finished = 0; + } + + t8_forest_set_user_data (bbox_forest, data); + + int refinement_pass = 0; + + t8_productionf ("starting refinement pass %i\n", refinement_pass); + + while (true) { + t8_forest_t forest_adapt; + t8_forest_init (&forest_adapt); + + t8_forest_set_adapt (forest_adapt, bbox_forest, t8_adapt_refine, 0); + t8_forest_ref (bbox_forest); + t8_forest_commit (forest_adapt); + + forest_data *new_data = forest_data_new (t8_forest_get_local_num_leaf_elements (forest_adapt)); + + t8_forest_set_user_data (forest_adapt, new_data); + t8_forest_iterate_replace (forest_adapt, bbox_forest, transfer_points); + + forest_data *old_data = static_cast (t8_forest_get_user_data (bbox_forest)); + + forest_data_destroy (old_data); + + t8_forest_unref (&bbox_forest); + + bbox_forest = forest_adapt; + + if (new_data->finished) { + break; + } + ++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_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); + // + // bbox_forest = adapted_forest; + // + // const t8_locidx_t leaf_count_after = t8_forest_get_local_num_leaf_elements (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; + // } + // + // t8_productionf ("adaptive refinement finished after %i pass(es)\n", refinement_pass + 1); + // if (!t8_forest_write_vtk (bbox_forest, "bounding_box_adapted")) { + // t8_productionf ("Could not write VTK file for forest"); + // }; + // + + forest_data *final_data = static_cast (t8_forest_get_user_data (bbox_forest)); + + t8_locidx_t new_tree_index = 0; + std::set mapped_old_tree_ids; + + const t8_locidx_t num_bbox_local_trees = t8_forest_get_num_local_trees (bbox_forest); + + for (t8_locidx_t bbox_itree = 0; bbox_itree < num_bbox_local_trees; ++bbox_itree) { + 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_locidx_t element_index = t8_forest_get_tree_element_offset (bbox_forest, bbox_itree) + ielement; + + const element_data &leaf_data = final_data->elem_data[element_index]; + + if (leaf_data.midpoints.empty ()) { + continue; + } + + if (leaf_data.midpoints.size () > 1) { + t8_productionf ("Warning: final bbox leaf still contains %lu tree centers. " + "Skipping this leaf.\n", + leaf_data.midpoints.size ()); + continue; + } + + const t8_gloidx_t old_tree_index = leaf_data.midpoints[0].first; + + if (!mapped_old_tree_ids.insert (old_tree_index).second) { + t8_productionf ("Warning: original tree %ld was already mapped. " + "Skipping duplicate occurrence.\n", + old_tree_index); + continue; + } + + tree_reindex[old_tree_index] = static_cast (new_tree_index); + + t8_productionf ("Original tree %u -> new SFC index %u\n", static_cast (old_tree_index), + static_cast (new_tree_index)); + + ++new_tree_index; + } + } + + 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_cmesh_trees)); + } + else { + t8_productionf ("successfully mapped all %u local tree(s)\n", static_cast (num_cmesh_trees)); + } + + forest_data_destroy (final_data); + + t8_forest_unref (&bbox_forest); + + t8_productionf ("tree reindexing finished\n"); + + return tree_reindex; +} + +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 = tree_reindex.at (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 = tree_reindex.at (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 = tree_reindex.at (old_id1); + + const t8_gloidx_t new_id2 = tree_reindex.at (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/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..b1d7103f1b --- /dev/null +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_tree_reindex.hxx @@ -0,0 +1,32 @@ +/* + 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, sc_MPI_Comm comm = sc_MPI_COMM_SELF); + +void +t8_cmesh_tree_perform_reindex_inplace (t8_stash_t &stash, const std::map &tree_reindex); diff --git a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_types.h b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_types.h index b7e47b1c41..952377b615 100644 --- a/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_types.h +++ b/src/t8_cmesh/t8_cmesh_internal/t8_cmesh_types.h @@ -99,6 +99,8 @@ typedef struct t8_cmesh int face_knowledge; /**< If partitioned the level of face knowledge that is expected. */ + int reindex_trees; /**< If nonzero the trees will be indexed according to their geometric locality */ + const t8_scheme_c *set_partition_scheme; /**< If the cmesh is to be partitioned according to a uniform level, the scheme that describes the refinement pattern. */ int8_t set_partition_level; /**< Non-negative if the cmesh should be partitioned from an already existing cmesh diff --git a/src/t8_config.h b/src/t8_config.h new file mode 100644 index 0000000000..1900bb4609 --- /dev/null +++ b/src/t8_config.h @@ -0,0 +1 @@ +//DO NOT COMMIT THIS \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f1e2f85152..65e99c9073 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -135,6 +135,8 @@ add_t8_cpp_test( NAME t8_gtest_cmesh_add_attributes_when_derive_parallel SOUR add_t8_cpp_test( NAME t8_gtest_cmesh_vertex_conn_tree_to_vertex_parallel SOURCES t8_cmesh/t8_gtest_cmesh_vertex_conn_tree_to_vertex.cxx ) add_t8_cpp_test( NAME t8_gtest_cmesh_vertex_conn_vertex_to_tree_parallel SOURCES t8_cmesh/t8_gtest_cmesh_vertex_conn_vertex_to_tree.cxx ) add_t8_cpp_test( NAME t8_gtest_cmesh_vertex_conn_serial SOURCES t8_cmesh/t8_gtest_cmesh_vertex_conn.cxx ) +add_t8_cpp_test( NAME t8_gtest_cmesh_tree_reindex_serial SOURCES t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx) +add_t8_cpp_test( NAME t8_gtest_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 ) add_t8_cpp_test( NAME t8_gtest_cmesh_bounding_box_serial SOURCES t8_cmesh/t8_gtest_cmesh_bounding_box.cxx ) diff --git a/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx b/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx new file mode 100644 index 0000000000..e53b0bff6c --- /dev/null +++ b/test/t8_cmesh/t8_gtest_cmesh_tree_reindex.cxx @@ -0,0 +1,290 @@ +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +struct t8_test_cmesh_tree_reindex: public testing::Test +{ + protected: + static constexpr std::size_t num_trees = 6; + + 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 + { + /* + * Same tetrahedral unit-cube decomposition as t8code's + * t8_cmesh_new_hypercube for T8_ECLASS_TET. + * + * Unit cube vertices: + * + * 0 = (0, 0, 0) + * 1 = (1, 0, 0) + * 2 = (0, 1, 0) + * 3 = (1, 1, 0) + * 4 = (0, 0, 1) + * 5 = (1, 0, 1) + * 6 = (0, 1, 1) + * 7 = (1, 1, 1) + */ + original_vertices = { /* + * Tree 0: vertices 0, 1, 5, 7 + */ + std::array { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0 }, + + /* + * Tree 1: vertices 0, 3, 1, 7 + */ + std::array { 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0 }, + + /* + * Tree 2: vertices 0, 2, 3, 7 + */ + std::array { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0 }, + + /* + * Tree 3: vertices 0, 6, 2, 7 + */ + std::array { 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0 }, + + /* + * Tree 4: vertices 0, 4, 6, 7 + */ + std::array { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0 }, + + /* + * Tree 5: vertices 0, 5, 4, 7 + */ + std::array { 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0 } + }; + + t8_cmesh_init (&cmesh); + ASSERT_NE (cmesh, nullptr); + + t8_cmesh_register_geometry (cmesh); + + for (t8_gloidx_t tree_id = 0; tree_id < static_cast (num_trees); ++tree_id) { + add_tet_tree (cmesh, tree_id, original_vertices[static_cast (tree_id)]); + } + + /* + * Same internal face joins as t8code's tetrahedral hypercube. + */ + t8_cmesh_set_join (cmesh, 0, 1, 2, 1, 0); + t8_cmesh_set_join (cmesh, 1, 2, 2, 1, 0); + t8_cmesh_set_join (cmesh, 2, 3, 2, 1, 0); + t8_cmesh_set_join (cmesh, 3, 4, 2, 1, 0); + t8_cmesh_set_join (cmesh, 4, 5, 2, 1, 0); + t8_cmesh_set_join (cmesh, 5, 0, 2, 1, 0); + } + + void + TearDown () override + { + if (cmesh != nullptr) { + t8_cmesh_unref (&cmesh); + cmesh = nullptr; + } + } + + t8_cmesh_t cmesh = nullptr; + std::array, num_trees> original_vertices; +}; + +TEST_F (t8_test_cmesh_tree_reindex, commit_reindexes_trees_successfully_and_correctly) +{ + t8_productionf ("Test started\n"); + + /* + * Build an identical uncommitted cmesh to compute the expected reindexing map. + * This cmesh is later committed without reindexing and written to VTK as a + * reference output. + */ + t8_cmesh_t expected_cmesh = nullptr; + t8_cmesh_init (&expected_cmesh); + + ASSERT_NE (expected_cmesh, nullptr); + ASSERT_FALSE (t8_cmesh_is_committed (expected_cmesh)); + + t8_cmesh_register_geometry (expected_cmesh); + + for (t8_gloidx_t tree_id = 0; tree_id < static_cast (num_trees); ++tree_id) { + add_tet_tree (expected_cmesh, tree_id, original_vertices[static_cast (tree_id)]); + } + + /* + * Same internal face joins as t8code's tetrahedral hypercube. + */ + t8_cmesh_set_join (expected_cmesh, 0, 1, 2, 1, 0); + t8_cmesh_set_join (expected_cmesh, 1, 2, 2, 1, 0); + t8_cmesh_set_join (expected_cmesh, 2, 3, 2, 1, 0); + t8_cmesh_set_join (expected_cmesh, 3, 4, 2, 1, 0); + t8_cmesh_set_join (expected_cmesh, 4, 5, 2, 1, 0); + t8_cmesh_set_join (expected_cmesh, 5, 0, 2, 1, 0); + + const std::map expected_reindex = t8_cmesh_reindex_tree (expected_cmesh, sc_MPI_COMM_SELF); + + ASSERT_FALSE (t8_cmesh_is_committed (expected_cmesh)); + + /* + * Check that the computed reindexing map is a valid bijection. + */ + ASSERT_EQ (expected_reindex.size (), num_trees); + + std::set old_tree_ids; + std::set new_tree_ids; + + bool reindex_is_identity = true; + + for (const auto &entry : expected_reindex) { + const t8_gloidx_t old_tree_id = entry.first; + const t8_gloidx_t new_tree_id = entry.second; + + EXPECT_GE (old_tree_id, 0); + EXPECT_LT (old_tree_id, static_cast (num_trees)); + + EXPECT_GE (new_tree_id, 0); + EXPECT_LT (new_tree_id, static_cast (num_trees)); + + old_tree_ids.insert (old_tree_id); + new_tree_ids.insert (new_tree_id); + + if (old_tree_id != new_tree_id) { + reindex_is_identity = false; + } + + t8_productionf ("Expected reindex: old global tree id %lli -> new global tree id %lli\n", + static_cast (old_tree_id), static_cast (new_tree_id)); + } + + EXPECT_EQ (old_tree_ids.size (), num_trees); + EXPECT_EQ (new_tree_ids.size (), num_trees); + + for (t8_gloidx_t tree_id = 0; tree_id < static_cast (num_trees); ++tree_id) { + EXPECT_EQ (old_tree_ids.count (tree_id), 1); + EXPECT_EQ (new_tree_ids.count (tree_id), 1); + } + + t8_productionf ("Reindexing is %s\n", reindex_is_identity ? "identity" : "non-identity"); + + /* + * Commit the reference cmesh without reindexing and write it to VTK. + */ + expected_cmesh->reindex_trees = 0; + + t8_cmesh_commit (expected_cmesh, sc_MPI_COMM_SELF); + + ASSERT_TRUE (t8_cmesh_is_committed (expected_cmesh)); + + t8_cmesh_vtk_write_file (expected_cmesh, "test_cmesh_tree_reindex_original"); + + /* + * Now commit the actual test cmesh with reindexing enabled. + */ + ASSERT_FALSE (t8_cmesh_is_committed (cmesh)); + ASSERT_NE (cmesh->stash, nullptr); + + cmesh->reindex_trees = 1; + + t8_cmesh_commit (cmesh, sc_MPI_COMM_SELF); + + ASSERT_TRUE (t8_cmesh_is_committed (cmesh)); + + t8_cmesh_vtk_write_file (cmesh, "test_cmesh_tree_reindex_reindexed"); + + for (const auto &entry : expected_reindex) { + const t8_gloidx_t old_tree_id = entry.first; + const t8_gloidx_t new_tree_id = entry.second; + + const t8_locidx_t new_local_tree_id = t8_cmesh_get_local_id (cmesh, new_tree_id); + + ASSERT_GE (new_local_tree_id, 0); + + EXPECT_EQ (t8_cmesh_get_global_id (cmesh, new_local_tree_id), new_tree_id); + + EXPECT_EQ (t8_cmesh_get_tree_class (cmesh, new_local_tree_id), T8_ECLASS_TET); + + double *actual_vertices = t8_cmesh_get_tree_vertices (cmesh, new_local_tree_id); + + ASSERT_NE (actual_vertices, nullptr); + + const std::array &expected_vertices = original_vertices[static_cast (old_tree_id)]; + + for (int icoord = 0; icoord < 12; ++icoord) { + EXPECT_DOUBLE_EQ (actual_vertices[icoord], expected_vertices[icoord]) + << "Mismatch for old tree id " << old_tree_id << ", new tree id " << new_tree_id << ", coordinate index " + << icoord; + } + + t8_productionf ("Verified old global tree id %lli -> new global tree id %lli\n", + static_cast (old_tree_id), static_cast (new_tree_id)); + } + + const std::array, 6> original_joins + = { std::array { 0, 1, 2, 1, 0 }, std::array { 1, 2, 2, 1, 0 }, + std::array { 2, 3, 2, 1, 0 }, std::array { 3, 4, 2, 1, 0 }, + std::array { 4, 5, 2, 1, 0 }, std::array { 5, 0, 2, 1, 0 } }; + + for (const auto &join : original_joins) { + const t8_gloidx_t old_tree_1 = static_cast (join[0]); + const t8_gloidx_t old_tree_2 = static_cast (join[1]); + + const int face_1 = join[2]; + const int face_2 = join[3]; + const int expected_orientation = join[4]; + + const t8_gloidx_t new_tree_1 = expected_reindex.at (old_tree_1); + const t8_gloidx_t new_tree_2 = expected_reindex.at (old_tree_2); + + const t8_locidx_t local_tree_1 = t8_cmesh_get_local_id (cmesh, new_tree_1); + + const t8_locidx_t local_tree_2 = t8_cmesh_get_local_id (cmesh, new_tree_2); + + ASSERT_GE (local_tree_1, 0); + ASSERT_GE (local_tree_2, 0); + + int dual_face = -1; + int orientation = -1; + + const t8_locidx_t neighbor_of_tree_1 + = t8_cmesh_get_face_neighbor (cmesh, local_tree_1, face_1, &dual_face, &orientation); + + ASSERT_GE (neighbor_of_tree_1, 0); + + EXPECT_EQ (t8_cmesh_get_global_id (cmesh, neighbor_of_tree_1), new_tree_2); + EXPECT_EQ (dual_face, face_2); + EXPECT_EQ (orientation, expected_orientation); + + dual_face = -1; + orientation = -1; + + const t8_locidx_t neighbor_of_tree_2 + = t8_cmesh_get_face_neighbor (cmesh, local_tree_2, face_2, &dual_face, &orientation); + + ASSERT_GE (neighbor_of_tree_2, 0); + + EXPECT_EQ (t8_cmesh_get_global_id (cmesh, neighbor_of_tree_2), new_tree_1); + EXPECT_EQ (dual_face, face_1); + + t8_productionf ("Verified reindexed join old trees %lli-%lli -> new trees %lli-%lli\n", + static_cast (old_tree_1), static_cast (old_tree_2), + static_cast (new_tree_1), static_cast (new_tree_2)); + } + + t8_cmesh_unref (&expected_cmesh); +}