From dbc50e1b0ed354800d3d6276b45cce77c2e92fd6 Mon Sep 17 00:00:00 2001 From: Alex Dutka <97711898+dutkalex@users.noreply.github.com> Date: Wed, 8 Oct 2025 11:22:12 +0200 Subject: [PATCH 1/8] fix empty struct warning --- src/t8_geometry/t8_geometry_hash.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/t8_geometry/t8_geometry_hash.hxx b/src/t8_geometry/t8_geometry_hash.hxx index 237b764df1..68084d87d8 100644 --- a/src/t8_geometry/t8_geometry_hash.hxx +++ b/src/t8_geometry/t8_geometry_hash.hxx @@ -32,13 +32,13 @@ #include #include -T8_EXTERN_C_BEGIN (); - /** Dummy tag for type trait usage of \ref t8_geometry_hash */ struct t8_geometry_hash_tag { }; +T8_EXTERN_C_BEGIN (); + /** Data type used for storing hash values of geometries. */ using t8_geometry_hash = T8Type; From 9e8ac7e6095ef5d2c6204937abaf9d8332196524 Mon Sep 17 00:00:00 2001 From: Thomas Spenke Date: Fri, 24 Oct 2025 15:53:35 +0200 Subject: [PATCH 2/8] Improvement: Move t8_cad.hxx/cxx to own folder src/t8_cad --- src/CMakeLists.txt | 4 ++-- src/{t8_data => t8_cad}/t8_cad.cxx | 2 +- src/{t8_data => t8_cad}/t8_cad.hxx | 0 .../t8_geometry_implementations/t8_geometry_cad.hxx | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) rename src/{t8_data => t8_cad}/t8_cad.cxx (99%) rename src/{t8_data => t8_cad}/t8_cad.hxx (100%) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d4d4e90ce2..d581a9a95c 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -108,12 +108,12 @@ if( T8CODE_ENABLE_OCC ) target_link_libraries( T8 PUBLIC ${OpenCASCADE_LIBRARIES} ) target_sources(T8 PRIVATE t8_geometry/t8_geometry_implementations/t8_geometry_cad.cxx - t8_data/t8_cad.cxx + t8_cad/t8_cad.cxx ) install( FILES t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx t8_geometry/t8_geometry_implementations/t8_geometry_cad.h - t8_data/t8_cad.hxx + t8_cad/t8_cad.hxx DESTINATION include ) endif() diff --git a/src/t8_data/t8_cad.cxx b/src/t8_cad/t8_cad.cxx similarity index 99% rename from src/t8_data/t8_cad.cxx rename to src/t8_cad/t8_cad.cxx index 61acf30452..f1a8f9d991 100644 --- a/src/t8_data/t8_cad.cxx +++ b/src/t8_cad/t8_cad.cxx @@ -21,7 +21,7 @@ */ #include -#include +#include #include #include #include diff --git a/src/t8_data/t8_cad.hxx b/src/t8_cad/t8_cad.hxx similarity index 100% rename from src/t8_data/t8_cad.hxx rename to src/t8_cad/t8_cad.hxx diff --git a/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx b/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx index ad87719632..03b6e0fc03 100644 --- a/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx +++ b/src/t8_geometry/t8_geometry_implementations/t8_geometry_cad.hxx @@ -33,7 +33,7 @@ #include #include #include -#include +#include #include #include #include From 56cf71f56e9abe69779b5b0e6d8abfadb566ad7f Mon Sep 17 00:00:00 2001 From: Johannes Holke Date: Wed, 12 Nov 2025 14:22:34 +0100 Subject: [PATCH 3/8] Add eclass string for invalid eclass --- src/t8_eclass.c | 2 +- src/t8_eclass.h | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/t8_eclass.c b/src/t8_eclass.c index 4d012d1624..257ac87f9e 100644 --- a/src/t8_eclass.c +++ b/src/t8_eclass.c @@ -65,7 +65,7 @@ const int t8_eclass_face_types[T8_ECLASS_COUNT][T8_ECLASS_MAX_FACES] = T8_ECLASS const int t8_eclass_boundary_count[T8_ECLASS_COUNT][T8_ECLASS_COUNT] = T8_ECLASS_BOUNDARY_COUNT_VALUES; -const char *t8_eclass_to_string[T8_ECLASS_COUNT] = T8_ECLASS_TO_STRING_VALUES; +const char *t8_eclass_to_string[T8_ECLASS_INVALID] = T8_ECLASS_TO_STRING_VALUES; int t8_eclass_count_boundary (t8_eclass_t theclass, int min_dim, int *per_eclass) diff --git a/src/t8_eclass.h b/src/t8_eclass.h index af9415da04..2e2728cac3 100644 --- a/src/t8_eclass.h +++ b/src/t8_eclass.h @@ -203,7 +203,7 @@ typedef enum t8_eclass { { 5, 8, 1, 4, 0, 0, 0, 0 } /* pyramid */ \ } -#define T8_ECLASS_TO_STRING_VALUES { "Vertex", "Line", "Quad", "Triangle", "Hex", "Tet", "Prism", "Pyramid" } +#define T8_ECLASS_TO_STRING_VALUES { "Vertex", "Line", "Quad", "Triangle", "Hex", "Tet", "Prism", "Pyramid", "Invalid" } /* clang-format on */ @@ -296,7 +296,7 @@ inline constexpr int t8_eclass_face_types[T8_ECLASS_COUNT][T8_ECLASS_MAX_FACES] inline constexpr int t8_eclass_boundary_count[T8_ECLASS_COUNT][T8_ECLASS_COUNT] = T8_ECLASS_BOUNDARY_COUNT_VALUES; /** For each eclass, the name of this class as a string */ -inline constexpr const char *t8_eclass_to_string[T8_ECLASS_COUNT] = T8_ECLASS_TO_STRING_VALUES; +inline constexpr const char *t8_eclass_to_string[T8_ECLASS_INVALID] = T8_ECLASS_TO_STRING_VALUES; } /* namespace t8cpp */ using namespace t8cpp; @@ -377,7 +377,7 @@ extern const int t8_eclass_face_types[T8_ECLASS_COUNT][T8_ECLASS_MAX_FACES]; extern const int t8_eclass_boundary_count[T8_ECLASS_COUNT][T8_ECLASS_COUNT]; /** For each eclass, the name of this class as a string */ -extern const char *t8_eclass_to_string[T8_ECLASS_COUNT]; +extern const char *t8_eclass_to_string[T8_ECLASS_INVALID]; #endif /* !__cplusplus */ From 208c91272a36eb90453c37f64032112cb0ebdebb Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer <49643115+sandro-elsweijer@users.noreply.github.com> Date: Mon, 17 Nov 2025 13:07:25 +0100 Subject: [PATCH 4/8] Update documentation link to Read the Docs --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 78f844477d..056cc69761 100644 --- a/README.md +++ b/README.md @@ -75,7 +75,7 @@ We provide a short guide to install t8code in our Wiki [Installation guide](http ### Documentation t8code uses [Doxygen](https://doxygen.nl/) to generate the code documentation. -You can find the documentation of our releases on the [t8code website](https://dlr-amr.github.io/t8code/pages/documentation.html). +You can find the documentation on [readthedocs](https://dlr-amr.github.io/t8code/pages/documentation.html). Follow the steps described in our Wiki [Documentation](https://github.com/DLR-AMR/t8code/wiki/Documentation) to create the documentation locally. From 33c8e148f9d4e8b6ff756f753a761cc70c1bf47a Mon Sep 17 00:00:00 2001 From: Johannes Holke Date: Tue, 28 Jan 2025 10:47:16 +0100 Subject: [PATCH 5/8] clean up custom assertion test header --- test/t8_gtest_custom_assertion.hxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/t8_gtest_custom_assertion.hxx b/test/t8_gtest_custom_assertion.hxx index b574bd6ede..fa5156852e 100644 --- a/test/t8_gtest_custom_assertion.hxx +++ b/test/t8_gtest_custom_assertion.hxx @@ -20,17 +20,17 @@ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ -/** \file t8_gtest_custom_assertion.cxx +/** \file t8_gtest_custom_assertion.hxx * Provide customized GoogleTest functions for improved error-output */ +#ifndef T8_GTEST_CUSTOM_ASSERTION_HXX +#define T8_GTEST_CUSTOM_ASSERTION_HXX + #include #include #include -#ifndef CUSTOM_ASSERTION_HXX -#define CUSTOM_ASSERTION_HXX - /** * \brief Test two elements for equality and print the elements if they aren't equal * @@ -104,4 +104,4 @@ vec_equality (const char *vec_1_expr, const char *vec_2_expr, const char *precis #define EXPECT_VEC_EQ(vec_1, vec_2, precision) EXPECT_PRED_FORMAT3 (vec_equality, (vec_1), (vec_2), (precision)) -#endif /* CUSTOM_ASSERTION_HXX */ +#endif /* T8_GTEST_CUSTOM_ASSERTION_HXX */ From a1a32bdfbdeab49545c0c5a08eabf977a3367481 Mon Sep 17 00:00:00 2001 From: Sandro Elsweijer <49643115+sandro-elsweijer@users.noreply.github.com> Date: Mon, 17 Nov 2025 14:54:58 +0100 Subject: [PATCH 6/8] Update documentation link to Read the Docs --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 056cc69761..0fdea22b35 100644 --- a/README.md +++ b/README.md @@ -75,7 +75,7 @@ We provide a short guide to install t8code in our Wiki [Installation guide](http ### Documentation t8code uses [Doxygen](https://doxygen.nl/) to generate the code documentation. -You can find the documentation on [readthedocs](https://dlr-amr.github.io/t8code/pages/documentation.html). +You can find the documentation on [readthedocs](https://t8code.readthedocs.io/en/latest/). Follow the steps described in our Wiki [Documentation](https://github.com/DLR-AMR/t8code/wiki/Documentation) to create the documentation locally. From 0edd339eb54e6c8b90e3d4ebe63f9c38768a8156 Mon Sep 17 00:00:00 2001 From: Johannes Holke Date: Tue, 18 Nov 2025 16:57:57 +0100 Subject: [PATCH 7/8] Add forest_neighbor folder with files for element and face neighbor --- src/t8_forest/t8_forest.cxx | 705 ------------------ src/t8_forest/t8_forest_general.h | 142 ---- .../t8_forest_element_face_neighbor.cxx | 215 ++++++ .../t8_forest_element_face_neighbor.h | 91 +++ .../t8_forest_leaf_face_neighbor.cxx | 558 ++++++++++++++ .../t8_forest_leaf_face_neighbor.h | 131 ++++ 6 files changed, 995 insertions(+), 847 deletions(-) create mode 100644 src/t8_forest/t8_forest_neighbor/t8_forest_element_face_neighbor.cxx create mode 100644 src/t8_forest/t8_forest_neighbor/t8_forest_element_face_neighbor.h create mode 100644 src/t8_forest/t8_forest_neighbor/t8_forest_leaf_face_neighbor.cxx create mode 100644 src/t8_forest/t8_forest_neighbor/t8_forest_leaf_face_neighbor.h diff --git a/src/t8_forest/t8_forest.cxx b/src/t8_forest/t8_forest.cxx index c841510a8c..79a4eafd1e 100644 --- a/src/t8_forest/t8_forest.cxx +++ b/src/t8_forest/t8_forest.cxx @@ -1417,711 +1417,6 @@ t8_forest_copy_trees (t8_forest_t forest, t8_forest_t from, int copy_elements) } } -/* TODO: This function seems to be untested. - * On Nov 7 2025 i found a bug in it that would have been caught by testing - * (the face number of the element was used as the cmesh tree face number, which - * is incorrect).*/ -t8_eclass_t -t8_forest_element_neighbor_eclass (const t8_forest_t forest, const t8_locidx_t ltreeid, const t8_element_t *elem, - const int face) -{ - /* Get a pointer to the tree to read its element class */ - const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, ltreeid); - const t8_scheme *scheme = t8_forest_get_scheme (forest); - if (!scheme->element_is_root_boundary (tree_class, elem, face)) { - /* The neighbor element is inside the current tree. */ - return tree_class; - } - /* We now compute the cmesh face neighbor class across the corresponding face. - * To do so, we first need to compute the face of the root tree corresponding to the - * face of the element. */ - const int tree_face = scheme->element_get_tree_face (tree_class, elem, face); - // Debug check if tree_face is a valid face of the tree. - // Note that if elem would not be a boundary element, the return value of - // element_get_tree_face could still be within these bounds, even though it does not make sense. - // We catch that case by checking element_is_root_boundary above. - T8_ASSERT (0 <= tree_face && tree_face < t8_eclass_num_faces[tree_class]); - - const t8_locidx_t cmesh_local_tree_id = t8_forest_ltreeid_to_cmesh_ltreeid (forest, ltreeid); - const t8_cmesh_t cmesh = t8_forest_get_cmesh (forest); - return t8_cmesh_get_tree_face_neighbor_eclass (cmesh, cmesh_local_tree_id, tree_face); -} - -/* TODO: If the forest has no ghosts, then skip the ghosts - parts. In that case, process boundary elements will have 0 neighbors. -*/ -t8_gloidx_t -t8_forest_element_face_neighbor (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *elem, t8_element_t *neigh, - const t8_eclass_t neigh_eclass, int face, int *neigh_face) -{ - /* Get a pointer to the tree to read its element class */ - const t8_eclass_t eclass = t8_forest_get_tree_class (forest, ltreeid); - const t8_scheme *scheme = t8_forest_get_scheme (forest); - if (neigh_eclass == eclass && scheme->element_get_face_neighbor_inside (eclass, elem, neigh, face, neigh_face)) { - /* The neighbor was constructed and is inside the current tree. */ - return t8_forest_global_tree_id (forest, ltreeid); - } - else { - /* The neighbor does not lie inside the current tree. The content of neigh is undefined right now. */ - t8_element_t *face_element; - int tree_neigh_face; - int orientation; - int is_smaller, eclass_compare; - - const t8_cmesh_t cmesh = forest->cmesh; - /* Get the scheme associated to the element class of the boundary element. */ - /* Compute the face of elem_tree at which the face connection is. */ - const int tree_face = scheme->element_get_tree_face (eclass, elem, face); - /* compute coarse tree id */ - const t8_locidx_t lctree_id = t8_forest_ltreeid_to_cmesh_ltreeid (forest, ltreeid); - T8_ASSERT (lctree_id >= 0); -#if T8_ENABLE_DEBUG - const bool cmesh_tree_is_local = t8_cmesh_treeid_is_local_tree (cmesh, lctree_id); - T8_ASSERT (cmesh_tree_is_local || t8_cmesh_treeid_is_ghost (cmesh, lctree_id)); -#endif - - const t8_locidx_t neighbor_ctreeid - = t8_cmesh_get_face_neighbor (cmesh, lctree_id, tree_face, &tree_neigh_face, &orientation); - - if (neighbor_ctreeid < 0) { - /* This face is a domain boundary. We do not need to continue */ - return -1; - } - - /* Get the eclass for the boundary */ - const t8_eclass_t boundary_class = (t8_eclass_t) t8_eclass_face_types[eclass][tree_face]; - /* Allocate the face element */ - scheme->element_new (boundary_class, 1, &face_element); - /* Compute the face element. */ - scheme->element_get_boundary_face (eclass, elem, face, face_element); - - /* We need to find out which face is the smaller one that is the one - * according to which the orientation was computed. - * face_a is smaller then face_b if either eclass_a < eclass_b - * or eclass_a = eclass_b and face_a < face_b. */ - /* -1 eclass < neigh_eclass, 0 eclass = neigh_eclass, 1 eclass > neigh_eclass */ - eclass_compare = t8_eclass_compare (eclass, neigh_eclass); - is_smaller = 0; - if (eclass_compare == -1) { - /* The face in the current tree is the smaller one */ - is_smaller = 1; - } - else if (eclass_compare == 1) { - /* The face in the other tree is the smaller one */ - is_smaller = 0; - } - else { - - T8_ASSERT (eclass_compare == 0); - /* Check if the face of the current tree has a smaller index then the face of the neighbor tree. */ - is_smaller = tree_face <= tree_neigh_face; - } - - /* We now transform the face element to the other tree. */ - const int sign - = t8_eclass_face_orientation[eclass][tree_face] == t8_eclass_face_orientation[neigh_eclass][tree_neigh_face]; - scheme->element_transform_face (boundary_class, face_element, face_element, orientation, sign, is_smaller); - /* And now we extrude the face to the new neighbor element */ - *neigh_face = scheme->element_extrude_face (neigh_eclass, face_element, neigh, tree_neigh_face); - /* Free the face_element */ - scheme->element_destroy (boundary_class, 1, &face_element); - - const t8_gloidx_t global_neigh_id = t8_cmesh_get_global_id (cmesh, neighbor_ctreeid); - - return global_neigh_id; - } -} - -t8_gloidx_t -t8_forest_element_half_face_neighbors (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *elem, - t8_element_t *neighs[], t8_eclass_t neigh_class, int face, int num_neighs, - int dual_faces[]) -{ - if (num_neighs <= 0) { - // There are no face neighbors. - // This case might happen and we need to catch it here - // before we use neighs[0] which might not be allocated. - return -1; - } - // Use the first allocated element temporarily as same level neighbors. - t8_element_t *same_level_neighbor = neighs[0]; - - // Compute the same level neighbor element. - int same_level_dual_face; - const t8_gloidx_t neighbor_tree = t8_forest_element_face_neighbor (forest, ltreeid, elem, same_level_neighbor, - neigh_class, face, &same_level_dual_face); - - if (neighbor_tree < 0) { - // No face neighbor exists - return -1; - } - // Double check that there is a neighbor tree. - // We expect the user to only call this function if neighbors exists (neigh_calls being an input argument). - T8_ASSERT (neighbor_tree >= 0 && neighbor_tree < t8_forest_get_num_global_trees (forest)); - - /* The scheme for the current forest */ - const t8_scheme *scheme = t8_forest_get_scheme (forest); - // Check that we are allowed to refine the neighbor element. - SC_CHECK_ABORT (scheme->element_get_level (neigh_class, same_level_neighbor) < t8_forest_get_maxlevel (forest), - "Trying to refine an element beyond its maximum allowed level."); - // Double check the number of neighbors. - T8_ASSERT (num_neighs - == scheme->element_get_num_face_children (neigh_class, same_level_neighbor, same_level_dual_face)); - - // Build the half face neighbors by constructing the children at the face. - scheme->element_get_children_at_face (neigh_class, same_level_neighbor, same_level_dual_face, neighs, num_neighs, - NULL); - - // We now need to compute the dual faces of the children. - // We do this with the scheme function - if (dual_faces != NULL) { - for (int iface_child = 0; iface_child < num_neighs; ++iface_child) { - dual_faces[iface_child] - = scheme->element_face_get_child_face (neigh_class, same_level_neighbor, same_level_dual_face, iface_child); - } - } - - return neighbor_tree; -} - -int -t8_forest_leaf_face_orientation (t8_forest_t forest, const t8_locidx_t ltreeid, const t8_scheme *scheme, - const t8_element_t *leaf, int face) -{ - int orientation = 0; - const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, ltreeid); - if (scheme->element_is_root_boundary (tree_class, leaf, face)) { - t8_cmesh_t cmesh = t8_forest_get_cmesh (forest); - t8_locidx_t ltreeid_in_cmesh = t8_forest_ltreeid_to_cmesh_ltreeid (forest, ltreeid); - int iface_in_tree = scheme->element_get_tree_face (tree_class, leaf, face); - t8_cmesh_get_face_neighbor (cmesh, ltreeid_in_cmesh, iface_in_tree, NULL, &orientation); - } - - return orientation; -} - -/** Internal data used for t8_forest_leaf_face_neighbors_iterate - * to buffer face neighbor information during leaf face neighbor search - * \ref t8_forest_leaf_face_neighbors_ext - * - * Given an element and a face, the search iterates through all leaves - * at that face and stores their information in this buffer. - * After the search the entries of the buffer are used and the - * search can start again with a clean buffer. -*/ -struct t8_lfn_user_data -{ - std::vector element_indices; /**< Element indices of the found neighbors. */ - std::vector dual_faces; /**< Dual faces of the found neighbors. */ - std::vector neighbors; /**< Pointers to the actual neighbor elements. */ -}; - -static int -t8_forest_leaf_face_neighbors_iterate (const t8_forest_t forest, const t8_locidx_t ltreeid, - [[maybe_unused]] const t8_element_t *element, const int face, const int is_leaf, - [[maybe_unused]] const t8_element_array_t *const leaf_elements, - const t8_locidx_t tree_leaf_index, void *user_data) -{ - // Output of iterate_faces: - // Array of indices in tree_leaves of all the face neighbor elements - // Assign pneighbor_leaves - // Assign dual_faces - // Assign pelement_indices - if (!is_leaf) { - // continue search until leaf level - return 1; - } - T8_ASSERT (is_leaf); - // Query whether this tree is a ghost and if so - // compute its id as a ghost tree ( 0 <= id < num_ghost_trees) - const bool is_ghost_tree = !t8_forest_tree_is_local (forest, ltreeid); - const t8_locidx_t adjusted_tree_id = !is_ghost_tree ? ltreeid : ltreeid - t8_forest_get_num_local_trees (forest); - T8_ASSERT (t8_forest_element_is_leaf_or_ghost (forest, element, adjusted_tree_id, is_ghost_tree)); - - struct t8_lfn_user_data *lfn_data = reinterpret_cast (user_data); - // face is the face of the considered leaf neighbor element and thus the - // corresponding dual face - t8_debugf ("Adding new face neighbor (leaf index %i) with dual face %i.\n", tree_leaf_index, face); - lfn_data->dual_faces.push_back (face); - // Compute the index of the element - const t8_locidx_t num_local_elements = t8_forest_get_local_num_leaf_elements (forest); - const t8_locidx_t tree_offset - = !is_ghost_tree ? t8_forest_get_tree_element_offset (forest, ltreeid) - : t8_forest_ghost_get_tree_element_offset (forest, adjusted_tree_id) + num_local_elements; - const t8_locidx_t element_index = tree_offset + tree_leaf_index; - lfn_data->element_indices.push_back (element_index); - // Add the pointer to the current element - const t8_element_t *&pnew_element = lfn_data->neighbors.emplace_back (); - if (!is_ghost_tree) { - pnew_element = t8_forest_get_leaf_element_in_tree (forest, ltreeid, tree_leaf_index); - } - else { - pnew_element = t8_forest_ghost_get_leaf_element (forest, adjusted_tree_id, tree_leaf_index); - } - return 1; -} - -/* - Set the proper return values of the leaf face neighbor computation - in case that no neighbors are found. - */ -static void -t8_forest_leaf_face_neighbors_set_no_neighbor_return_value (const t8_element_t **pneighbor_leaves[], int *dual_faces[], - int *num_neighbors, t8_locidx_t **pelement_indices, - t8_gloidx_t *gneigh_tree) -{ - *dual_faces = NULL; - *num_neighbors = 0; - *pelement_indices = NULL; - if (pneighbor_leaves) { - *pneighbor_leaves = NULL; - } - if (gneigh_tree != NULL) { - *gneigh_tree = -1; - } -} - -/* TODO: If the forest has no ghosts, then skip the ghosts - parts. In that case, process boundary elements will have 0 neighbors. -*/ -void -t8_forest_leaf_face_neighbors_ext (const t8_forest_t forest, const t8_locidx_t ltreeid, - const t8_element_t *leaf_or_ghost, const t8_element_t **pneighbor_leaves[], - const int face, int *dual_faces[], int *num_neighbors, - t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass, t8_gloidx_t *gneigh_tree, - int *orientation) -{ - /* We compute all face neighbor leaf elements of E via the following strategy: - * - Compute the same level face neighbor N - * - The neighbor tree could be a local tree or ghost (or both), - * for each variant get the leaf array of the neighbor tree and search in it: - * - Search for the first leaf element L overlapping with N. - * If it exists, it is either an ancestor or descendant of N (or N itself which is included in both definitions). - * If it does not exist, there are not leaf face neighbors in this tree. - * - If L is an ancestor of N (i.e. Level(L) <= Level(N)) then L is the only face neighbor. - * - Otherwise (Level(L) > Level (N)) we use a recursive face search across N's neighbor face, - * adding all leaf elements on the face to the face neighbors. - * This search will require L (more precise its position in the tree leaf array) as input. - **/ - - T8_ASSERT (t8_forest_is_committed (forest)); - T8_ASSERT (pelement_indices != NULL); - T8_ASSERT (dual_faces != NULL); - T8_ASSERT (num_neighbors != NULL); - -#if T8_ENABLE_DEBUG - const bool tree_is_local = t8_forest_tree_is_local (forest, ltreeid); - if (tree_is_local) { - T8_ASSERT (t8_forest_element_is_leaf (forest, leaf_or_ghost, ltreeid)); - } - else { - const t8_locidx_t local_ghost_treeid = ltreeid - t8_forest_get_num_local_trees (forest); - T8_ASSERT (t8_forest_element_is_ghost (forest, leaf_or_ghost, local_ghost_treeid)); - } -#endif - SC_CHECK_ABORT (!forest->incomplete_trees, "Leaf face neighbor is not supported for " - "forests with deleted elements.\n"); - - const t8_scheme *scheme = t8_forest_get_scheme (forest); - - if (orientation) { - // Compute the orientation of the face neighbor connection - *orientation = t8_forest_leaf_face_orientation (forest, ltreeid, scheme, leaf_or_ghost, face); - } - - /* At first we compute the same lave face neighbor element of leaf. For this, we need the - * neighbor tree's eclass and scheme. */ - const t8_eclass_t neigh_class = t8_forest_element_neighbor_eclass (forest, ltreeid, leaf_or_ghost, face); - if (neigh_class == T8_ECLASS_INVALID) { - // No face neighbor exists. - t8_forest_leaf_face_neighbors_set_no_neighbor_return_value (pneighbor_leaves, dual_faces, num_neighbors, - pelement_indices, gneigh_tree); - return; - } - if (pneigh_eclass != NULL) { - *pneigh_eclass = neigh_class; - } - - // Compute the same level face neighbor - t8_element_t *same_level_neighbor; - scheme->element_new (neigh_class, 1, &same_level_neighbor); - int same_level_neighbor_dual_face; - const t8_gloidx_t computed_gneigh_tree = t8_forest_element_face_neighbor ( - forest, ltreeid, leaf_or_ghost, same_level_neighbor, neigh_class, face, &same_level_neighbor_dual_face); - t8_debugf ("Computed same level neighbor with dual face %i\n", same_level_neighbor_dual_face); - - if (computed_gneigh_tree < 0) { - // There is no face neighbor across this face - scheme->element_destroy (neigh_class, 1, &same_level_neighbor); - t8_forest_leaf_face_neighbors_set_no_neighbor_return_value (pneighbor_leaves, dual_faces, num_neighbors, - pelement_indices, gneigh_tree); - return; - } - - // The neighbor leaves could be distributed across a local tree and a ghost - // tree. We thus possibly need to search in two different arrays. - // We store these in a vector and iterate over the entries. - // The leaf arrays themself do not store any information about their tree, - // whether it is local or ghost. - // We thus need to add this info and hence store a pair of element array and - // a bool that is true if and only if the element array corresponds to a ghost tree. - using neighbor_leaf_array = std::pair; - - // We compute the owners of the first and last face descendant. - // If the current rank is in between then the local process might have neighbor elements - // and we search the local tree. - // If other processes are in the interval of owners (lower_bound < q != p < upper_bound), - // then we (additionally or alone) search the ghost tree. - int face_owners_lower_bound = 0; - int face_owners_upper_bound = forest->mpisize - 1; - const int mpirank = forest->mpirank; - t8_forest_element_owners_at_face_bounds (forest, computed_gneigh_tree, same_level_neighbor, neigh_class, - same_level_neighbor_dual_face, &face_owners_lower_bound, - &face_owners_upper_bound); - - std::vector leaf_arrays; - - const t8_locidx_t local_neighbor_tree = t8_forest_get_local_id (forest, computed_gneigh_tree); - if (face_owners_lower_bound <= mpirank && mpirank <= face_owners_upper_bound) { - // Add the local neighbor tree's elements to the search array. - // Compute the local id of the neighbor tree and check if it is a local tree - t8_debugf ("Adding local tree to search.\n"); - if (0 <= local_neighbor_tree) { - // The neighbor tree is a local tree and hence there may be local neighbor elements. - const t8_element_array_t *tree_leaves = t8_forest_tree_get_leaf_elements (forest, local_neighbor_tree); - if (tree_leaves != nullptr) { - neighbor_leaf_array *leaf_array = new neighbor_leaf_array (tree_leaves, false); - leaf_arrays.push_back (leaf_array); - } - } - } - - if (forest->ghosts != NULL) { - if (face_owners_lower_bound != mpirank || face_owners_upper_bound != mpirank) { - // Add the neighbor tree ghost elements to the search array - t8_debugf ("Adding ghost tree to search.\n"); - const t8_locidx_t local_neighbor_ghost_treeid = t8_forest_ghost_get_ghost_treeid (forest, computed_gneigh_tree); - if (local_neighbor_ghost_treeid >= 0) { - // The neighbor tree is also a ghost tree and face neighbors of our element might - // be ghost elements. - // We add the ghost elements of that tree to our search array. - const t8_element_array_t *ghost_leaves - = t8_forest_ghost_get_tree_leaf_elements (forest, local_neighbor_ghost_treeid); - if (ghost_leaves != nullptr) { - neighbor_leaf_array *leaf_array = new neighbor_leaf_array (ghost_leaves, true); - leaf_arrays.push_back (leaf_array); - } - } - } - } - - struct t8_lfn_user_data user_data; - - // Now we iterate over the leaf arrays of the neighbor tree - // or neighbor ghost tree and find all leaf face neighbors of the element. - *num_neighbors = 0; - // Since we use REALLOC later to allocate memory of the following - // three pointers, we have to set them to NULL manually. - // This will trigger REALLOC to allocate the memory in the initial call. - // Not setting them to NULL but keeping them possibly uninitialized, will - // call REALLOC on uninitialized memory and result in memory errors. - if (pneighbor_leaves != NULL) { - /* Only set *pneighbor_leaves if a computation is desired. */ - *pneighbor_leaves = NULL; - } - - *pelement_indices = NULL; - *dual_faces = NULL; - for (auto &leaf_array : leaf_arrays) { - auto &tree_leaves = leaf_array->first; - const bool leaf_array_is_ghost = leaf_array->second; - T8_ASSERT (tree_leaves != NULL); - const t8_element_t *first_descendant; - /* - * Compute the index of the first leaf in tree_leaves that is an ancestor or descendant of - * the same_level_neighbor (might be the neighbor itself). - * Such an element might not exist in which case there are no neighbors in this tree_leaves - * array. - */ - const t8_locidx_t first_leaf_index - = t8_forest_bin_search_first_descendant_ancestor (tree_leaves, same_level_neighbor, &first_descendant); - - if (first_leaf_index >= 0) { - T8_ASSERT (first_descendant != nullptr); - const int neighbor_level = scheme->element_get_level (neigh_class, same_level_neighbor); - const int first_desc_level = scheme->element_get_level (neigh_class, first_descendant); - /* If the same level neighbor is coarser than the first found leaf, then - * we iterate over the faces of the same level neighbor. - * Otherwise, there is only one face neighbor, the first_descendant. - * We will do the iteration over the first_descendant nevertheless, but it will stop immediately. - */ - T8_ASSERT (neighbor_level >= 0); - T8_ASSERT (first_desc_level >= 0); - const bool neighbor_unique = first_desc_level <= neighbor_level; - const t8_element_t *search_this_element = neighbor_unique ? first_descendant : same_level_neighbor; - t8_debugf ("[H] Starting face search. neigh level %i, desc level %i\n", neighbor_level, first_desc_level); - - int temp_dual_face; - if (neighbor_unique && first_desc_level <= neighbor_level) { - temp_dual_face = scheme->element_face_get_ancestor_face (neigh_class, same_level_neighbor, first_desc_level, - same_level_neighbor_dual_face); - } // end if neighbor_unique - - const int search_element_dual_face = neighbor_unique ? temp_dual_face : same_level_neighbor_dual_face; - t8_debugf ("\tSearch dual face is %i\n", search_element_dual_face); - - // There may be face neighbors in this leaf array. - - // We need to restrict the array such that it contains only elements inside the search element. - // Thus, we create a new view containing all elements starting at first_leaf_index. - t8_element_array_t reduced_leaves; - if (neighbor_unique) { - t8_debugf ("Starting search with element indices %i to %i (including).\n", first_leaf_index, first_leaf_index); - t8_element_array_init_view (&reduced_leaves, tree_leaves, first_leaf_index, 1); - } - else { - /* We need to compute the first element that is not longer contained in the same_level_neighbor. - * To do so, we compute the successor of the same_level_neighbor and do - * an upper search for it in the leaf array. - * The found element (if existing) is the first leaf that is not a descendant of same_level_neighbor. */ - /* The successor might not exist because the same level neighbor is the last - * element of its level in the tree. - * To identify this case, we check whether the last leaf of the tree is an - * ancestor of same_level_neighbor. If so, then it is automatically the last leaf - * that we need to check. - * If not, we build the successor of same_level_neighbor. */ - const t8_locidx_t leaf_count = t8_element_array_get_count (tree_leaves); - const t8_element_t *last_leaf = t8_element_array_index_locidx (tree_leaves, leaf_count - 1); - T8_ASSERT (last_leaf != NULL); - t8_locidx_t last_search_element_index = -1; - if (scheme->element_is_ancestor (neigh_class, same_level_neighbor, last_leaf)) { - last_search_element_index = leaf_count - 1; - } - else { - - t8_element_t *successor; - scheme->element_new (neigh_class, 1, &successor); - scheme->element_construct_successor (neigh_class, same_level_neighbor, successor); - const int successor_level = scheme->element_get_level (neigh_class, successor); - const t8_linearidx_t successor_id = scheme->element_get_linear_id (neigh_class, successor, successor_level); - scheme->element_destroy (neigh_class, 1, &successor); - const t8_locidx_t upper_search_for_successor - = t8_forest_bin_search_upper (tree_leaves, successor_id, successor_level); - // The first index of a non descendant is the found element or the end of the array - // if no element was found. - // The last index in our search range is 1 less. - last_search_element_index = upper_search_for_successor < 0 ? leaf_count - 1 : upper_search_for_successor - 1; - } - const size_t reduced_leaf_count = last_search_element_index - first_leaf_index + 1; - T8_ASSERT (reduced_leaf_count > 0); - t8_debugf ("Starting search with element indices %i to %i (including).\n", first_leaf_index, - last_search_element_index); - t8_element_array_init_view (&reduced_leaves, tree_leaves, first_leaf_index, reduced_leaf_count); - } - // Iterate over all leaves at the face and collect them as neighbors. - const t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (forest); - // Compute the local or ghost tree id depending on whether this leaf array corresponds to a local - // tree or ghost tree. - const t8_locidx_t face_iterate_tree_id - = leaf_array_is_ghost ? t8_forest_ghost_get_ghost_treeid (forest, computed_gneigh_tree) + num_local_trees - : local_neighbor_tree; - t8_forest_iterate_faces (forest, face_iterate_tree_id, search_this_element, search_element_dual_face, - &reduced_leaves, first_leaf_index, t8_forest_leaf_face_neighbors_iterate, &user_data); - // Output of iterate_faces: - // Array of indices in tree_leaves of all the face neighbor elements - // Assign pneighbor_leaves - // Assign dual_faces - // Assign pelement_indices - // (all as growing std::vectors, resp t8_element_array) - - // num_neighbors counts the already inserted neighbors before this tree - // num_neighbors_current_tree counts the neighbors added in this tree - // total_num_neighbors temporarily counts all inserted neighbors, including this tree - const int num_neighbors_current_tree = user_data.neighbors.size (); - const int total_num_neighbors = *num_neighbors + num_neighbors_current_tree; - t8_debugf ("Found %i neighbors in tree. Adding up to %i total neighbors.\n", num_neighbors_current_tree, - total_num_neighbors); - // Copy neighbor element pointers - if (num_neighbors_current_tree > 0) { - if (pneighbor_leaves != NULL) { - *pneighbor_leaves = T8_REALLOC (*pneighbor_leaves, const t8_element_t *, total_num_neighbors); - T8_ASSERT (*pneighbor_leaves != NULL); - // Copy the pointers to pneighbor_leaves - for (t8_locidx_t ielem = 0; ielem < num_neighbors_current_tree; ++ielem) { - (*pneighbor_leaves)[ielem] = user_data.neighbors.data ()[*num_neighbors + ielem]; - } - } - // Copy element indices - *pelement_indices = T8_REALLOC (*pelement_indices, t8_locidx_t, total_num_neighbors); - T8_ASSERT (*pelement_indices != NULL); - memcpy (*pelement_indices + *num_neighbors, user_data.element_indices.data () + *num_neighbors, - num_neighbors_current_tree * sizeof (t8_locidx_t)); - // Copy dual face - *dual_faces = T8_REALLOC (*dual_faces, int, total_num_neighbors); - T8_ASSERT (*dual_faces != NULL); - memcpy (*dual_faces + *num_neighbors, user_data.dual_faces.data () + *num_neighbors, - num_neighbors_current_tree * sizeof (int)); - *num_neighbors = total_num_neighbors; - } - } // End if neighbors exist (first_leaf_index > 0) - // clean up memory allocated with new - delete leaf_array; - } - scheme->element_destroy (neigh_class, 1, &same_level_neighbor); -#if T8_ENABLE_DEBUG - // Debugging checks - if (tree_is_local && forest->ghosts != NULL) { - // For local elements we must have found face neighbors by now. - T8_ASSERT (*num_neighbors > 0); - } - // All neighbor elements must be valid - if (pneighbor_leaves != NULL) { - for (int ineigh = 0; ineigh < *num_neighbors; ++ineigh) { - T8_ASSERT (scheme->element_is_valid (neigh_class, (*pneighbor_leaves)[ineigh])); - t8_debugf ("Face neighbor %p is valid.\n", (void *) (*pneighbor_leaves)[ineigh]); - scheme->element_debug_print (neigh_class, (*pneighbor_leaves)[ineigh]); - } - } -#endif // T8_ENABLE_DEBUG - // If no neighbors are found, set the proper return values. - if (*num_neighbors == 0) { - t8_debugf ("Found no neighbors\n"); - t8_forest_leaf_face_neighbors_set_no_neighbor_return_value (pneighbor_leaves, dual_faces, num_neighbors, - pelement_indices, gneigh_tree); - T8_ASSERT (*dual_faces == NULL); - T8_ASSERT (*pelement_indices == NULL); - T8_ASSERT (pneighbor_leaves == NULL || *pneighbor_leaves == NULL); - T8_ASSERT (gneigh_tree == NULL || *gneigh_tree == -1); - } - - if (gneigh_tree != NULL) { - *gneigh_tree = computed_gneigh_tree; - } -} - -void -t8_forest_leaf_face_neighbors (const t8_forest_t forest, const t8_locidx_t ltreeid, const t8_element_t *leaf, - const t8_element_t **pneighbor_leaves[], const int face, int *dual_faces[], - int *num_neighbors, t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass) -{ - t8_forest_leaf_face_neighbors_ext (forest, ltreeid, leaf, pneighbor_leaves, face, dual_faces, num_neighbors, - pelement_indices, pneigh_eclass, NULL, NULL); -} - -t8_locidx_t -t8_forest_same_level_leaf_face_neighbor_index (const t8_forest_t forest, const t8_locidx_t element_index, - const int face_index, const t8_gloidx_t global_treeid, int *dual_face) -{ - const t8_locidx_t num_local_elements = t8_forest_get_local_num_leaf_elements (forest); -#if T8_ENABLE_DEBUG - const t8_locidx_t num_ghosts = t8_forest_get_num_ghosts (forest); - T8_ASSERT (0 <= element_index && element_index < num_local_elements + num_ghosts); -#endif - const bool is_local = element_index < num_local_elements; - - t8_locidx_t local_tree; - t8_locidx_t element_index_in_tree; - const t8_element_t *element; - if (is_local) { - local_tree = t8_forest_get_local_id (forest, global_treeid); - element_index_in_tree = element_index - t8_forest_get_tree_element_offset (forest, local_tree); - element = t8_forest_get_leaf_element_in_tree (forest, local_tree, element_index_in_tree); - } - else { - local_tree = t8_forest_ghost_get_ghost_treeid (forest, global_treeid); - const t8_locidx_t ghost_offset_in_tree = t8_forest_ghost_get_tree_element_offset (forest, local_tree); - element_index_in_tree = element_index - num_local_elements - ghost_offset_in_tree; - element = t8_forest_ghost_get_leaf_element (forest, local_tree, element_index_in_tree); - local_tree += t8_forest_get_num_local_trees (forest); - } - - int *dual_faces; - int num_neighbors = 0; - t8_locidx_t *element_indices; - t8_eclass_t neigh_class; - - t8_debugf ("Same level leaf neighbor for index %i. Which is %s element %i in tree %i.\n", element_index, - element_index < num_local_elements ? "local" : "ghost", element_index_in_tree, local_tree); - - t8_forest_leaf_face_neighbors (forest, local_tree, element, NULL, face_index, &dual_faces, &num_neighbors, - &element_indices, &neigh_class); - - T8_ASSERT (num_neighbors == 0 || num_neighbors == 1); - - if (num_neighbors == 0) { - *dual_face = -1; - return -1; - } - - *dual_face = dual_faces[0]; - const t8_locidx_t neigh_index = element_indices[0]; - - T8_FREE (element_indices); - T8_FREE (dual_faces); - - return neigh_index; -} - -void -t8_forest_print_all_leaf_neighbors (t8_forest_t forest) -{ - t8_locidx_t ltree, ielem; - const t8_element_t **neighbor_leaves; - int iface, num_neighbors, ineigh; - t8_eclass_t eclass, neigh_eclass; - const t8_scheme *scheme = t8_forest_get_scheme (forest); - t8_locidx_t *element_indices; - int *dual_faces; - char buffer[BUFSIZ]; - int allocate_first_desc = 0, allocate_tree_offset = 0; - int allocate_el_offset = 0; - - if (forest->tree_offsets == NULL) { - allocate_tree_offset = 1; - t8_forest_partition_create_tree_offsets (forest); - } - if (forest->global_first_desc == NULL) { - allocate_first_desc = 1; - t8_forest_partition_create_first_desc (forest); - } - if (forest->element_offsets == NULL) { - allocate_el_offset = 1; - t8_forest_partition_create_offsets (forest); - } - for (ielem = 0; ielem < t8_forest_get_local_num_leaf_elements (forest); ielem++) { - /* Get a pointer to the ielem-th element, its eclass, treeid and scheme */ - const t8_element_t *leaf = t8_forest_get_leaf_element (forest, ielem, <ree); - eclass = t8_forest_get_tree_class (forest, ltree); - /* Iterate over all faces */ - for (iface = 0; iface < scheme->element_get_num_faces (eclass, leaf); iface++) { - t8_forest_leaf_face_neighbors (forest, ltree, leaf, &neighbor_leaves, iface, &dual_faces, &num_neighbors, - &element_indices, &neigh_eclass); - t8_debugf ("Element %li across face %i has %i leaf neighbors (with dual faces).\n", (long) ielem, iface, - num_neighbors); - snprintf (buffer, BUFSIZ, "\tIndices:\t"); - for (ineigh = 0; ineigh < num_neighbors; ineigh++) { - snprintf (buffer + strlen (buffer), BUFSIZ - strlen (buffer), "%li (%i) ", (long) element_indices[ineigh], - dual_faces[iface]); - } - t8_debugf ("%s\n", buffer); - if (num_neighbors > 0) { - T8_FREE (element_indices); - T8_FREE (neighbor_leaves); - T8_FREE (dual_faces); - } - } - } - if (allocate_tree_offset) { - t8_shmem_array_destroy (&forest->tree_offsets); - } - if (allocate_first_desc) { - t8_shmem_array_destroy (&forest->global_first_desc); - } - if (allocate_el_offset) { - t8_shmem_array_destroy (&forest->element_offsets); - } -} - int t8_forest_tree_is_local (const t8_forest_t forest, const t8_locidx_t local_tree) { diff --git a/src/t8_forest/t8_forest_general.h b/src/t8_forest/t8_forest_general.h index 09d45157b3..f4b3716275 100644 --- a/src/t8_forest/t8_forest_general.h +++ b/src/t8_forest/t8_forest_general.h @@ -550,112 +550,6 @@ int t8_forest_element_is_leaf_or_ghost (const t8_forest_t forest, const t8_element_t *element, const t8_locidx_t local_tree, const int check_ghost); -/** Compute the leaf face orientation at given face in a forest. - * \param [in] forest The forest. Must have a valid ghost layer. - * \param [in] ltreeid A local tree id. - * \param [in] scheme The eclass scheme of the element. - * \param [in] leaf A leaf in tree \a ltreeid of \a forest. - * \param [in] face The index of the face across which the face neighbors - * are searched. - * \return Face orientation encoded as integer. - * - * For more information about the encoding of face orientation refer to \ref t8_cmesh_get_face_neighbor. - */ -int -t8_forest_leaf_face_orientation (t8_forest_t forest, const t8_locidx_t ltreeid, const t8_scheme_c *scheme, - const t8_element_t *leaf, const int face); - -/** Compute the leaf face neighbors of a forest leaf element or ghost leaf. - * \param [in] forest The forest. - * \param [in] ltreeid A local tree id (could also be a ghost tree). 0 <= \a ltreeid < num_local trees+num_ghost_trees - * \param [in] leaf A leaf in tree \a ltreeid of \a forest. - * \param [out] pneighbor_leaves Unallocated on input. On output the neighbor - * leaves are stored here. - * \param [in] face The index of the face across which the face neighbors - * are searched. - * \param [out] dual_faces On output the face id's of the neighboring elements' faces. - * \param [out] num_neighbors On output the number of neighbor leaves. - * \param [out] pelement_indices Unallocated on input. On output the element indices - * of the neighbor leaves are stored here. - * 0, 1, ... num_local_el - 1 for local leaves and - * num_local_el , ... , num_local_el + num_ghosts - 1 for ghosts. - * \param [out] pneigh_eclass On output the eclass of the neighbor elements. - * \note If there are no face neighbors, then *pneighbor_leaves = NULL, num_neighbors = 0, - * and *pelement_indices = NULL on output. - * \note \a forest must be committed before calling this function. - * \note If \a forest does not have a ghost layer then leaf elements at the process boundaries have 0 neighbors. (The function output for leaf elements then depends on the parallel partition.) - * \note Important! This routine allocates memory which must be freed. Do it like this: - * - * if (num_neighbors > 0) { - * T8_FREE (pneighbor_leaves); - * T8_FREE (pelement_indices); - * T8_FREE (dual_faces); - * } - * - */ -void -t8_forest_leaf_face_neighbors (const t8_forest_t forest, const t8_locidx_t ltreeid, const t8_element_t *leaf, - const t8_element_t **pneighbor_leaves[], const int face, int *dual_faces[], - int *num_neighbors, t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass); - -/** Like \ref t8_forest_leaf_face_neighbors but also provides information about the global neighbors and the orientation. - * \param [in] forest The forest. Must have a valid ghost layer. - * \param [in] ltreeid A local tree id (could also be a ghost tree). 0 <= \a ltreeid < num_local trees+num_ghost_trees - * \param [in] leaf A leaf in tree \a ltreeid of \a forest. - * \param [out] pneighbor_leaves Unallocated on input. On output the neighbor - * leaves are stored here. - * \param [in] face The index of the face across which the face neighbors - * are searched. - * \param [out] dual_faces On output the face id's of the neighboring elements' faces. - * \param [out] num_neighbors On output the number of neighbor leaves. - * \param [out] pelement_indices Unallocated on input. On output the element indices - * of the neighbor leaves are stored here. - * 0, 1, ... num_local_el - 1 for local leaves and - * num_local_el , ... , num_local_el + num_ghosts - 1 for ghosts. - * \param [out] pneigh_eclass On output the eclass of the neighbor elements. - * \param [out] gneigh_tree The global tree IDs of the neighbor trees. - * \param [out] orientation If not NULL on input, the face orientation is computed and stored here. - * Thus, if the face connection is an inter-tree connection the orientation of the tree-to-tree connection is stored. - * Otherwise, the value 0 is stored. - * All other parameters and behavior are identical to \ref t8_forest_leaf_face_neighbors. - * \note If there are no face neighbors, then *pneighbor_leaves = NULL, num_neighbors = 0, - * and *pelement_indices = NULL on output. - * \note \a forest must be committed before calling this function. - * - * \note Important! This routine allocates memory which must be freed. Do it like this: - * - * if (num_neighbors > 0) { - * T8_FREE (pneighbor_leaves); - * T8_FREE (pelement_indices); - * T8_FREE (dual_faces); - * } - * - */ -void -t8_forest_leaf_face_neighbors_ext (const t8_forest_t forest, const t8_locidx_t ltreeid, - const t8_element_t *leaf_or_ghost, const t8_element_t **pneighbor_leaves[], - const int face, int *dual_faces[], int *num_neighbors, - t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass, t8_gloidx_t *gneigh_tree, - int *orientation); - -/** Given a leaf element or ghost index in "all local elements + ghosts" enumeration - * compute the index of the face neighbor of the element - provided that only one or no - * face neighbors exists. - * HANDLE WITH CARE. DO NOT CALL IF THE FOREST IS ADAPTED. - * - * \param[in] forest The forest. Must be committed. - * \param[in] element_index Index of an element in \a forest. Must have only one or no facen neighbors across the given face. - * 0 <= \a element_index < num_local_elements + num_ghosts - * \param[in] face_index Index of a face of \a element. - * \param[in] global_treeid Global index of the tree that contains \a element. - * \param[out] dual_face Return value, the dual_face index of the face neighbor. - * \return The index of the face neighbor leaf (local element or ghost). - * \note Do not call if you are unsure about the number of face neighbors. In particular if the forest is adapted and not uniform. - */ -t8_locidx_t -t8_forest_same_level_leaf_face_neighbor_index (const t8_forest_t forest, const t8_locidx_t element_index, - const int face_index, const t8_gloidx_t global_treeid, int *dual_face); - /** Exchange ghost information of user defined element data. * \param [in] forest The forest. Must be committed. * \param [in] element_data An array of length num_local_elements + num_ghosts @@ -844,42 +738,6 @@ t8_forest_get_first_local_leaf_element_id (t8_forest_t forest); const t8_scheme_c * t8_forest_get_scheme (const t8_forest_t forest); -/** Return the eclass of the tree in which a face neighbor of a given element or ghost - * lies. - * \param [in] forest A committed forest. - * \param [in] ltreeid The local tree or ghost tree in which the element lies. 0 <= \a ltreeid < num_local_trees + num_ghost_trees - * \param [in] elem An element or ghost in the tree \a ltreeid. - * \param [in] face A face number of \a elem. - * \return The eclass of the local tree or ghost tree that - * is face neighbor of \a elem across \a face. - * T8_ECLASS_INVALID if no neighbor exists. - */ -t8_eclass_t -t8_forest_element_neighbor_eclass (const t8_forest_t forest, const t8_locidx_t ltreeid, const t8_element_t *elem, - const int face); - -/** Construct the face neighbor of an element, possibly across tree boundaries. - * Returns the global tree-id of the tree in which the neighbor element lies in. - * - * \param [in] forest The forest. - * \param [in] ltreeid The local tree in which the element lies. - * \param [in] elem The element to be considered. - * \param [in,out] neigh On input an allocated element of the scheme of the - * face_neighbors eclass. - * On output, this element's data is filled with the - * data of the face neighbor. If the neighbor does not exist - * the data could be modified arbitrarily. - * \param [in] neigh_eclass The eclass of \a neigh. - * \param [in] face The number of the face along which the neighbor should be constructed. - * \param [out] neigh_face The number of the face viewed from perspective of \a neigh. - * \return The global tree-id of the tree in which \a neigh is in. - * -1 if there exists no neighbor across that face. Domain boundary. - * -2 if the neighbor is not in a local tree or ghost tree. Process/Ghost boundary. - */ -t8_gloidx_t -t8_forest_element_face_neighbor (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *elem, t8_element_t *neigh, - const t8_eclass_t neigh_eclass, int face, int *neigh_face); - /** * TODO: Can be removed since it is unused. * diff --git a/src/t8_forest/t8_forest_neighbor/t8_forest_element_face_neighbor.cxx b/src/t8_forest/t8_forest_neighbor/t8_forest_element_face_neighbor.cxx new file mode 100644 index 0000000000..0ce32776b4 --- /dev/null +++ b/src/t8_forest/t8_forest_neighbor/t8_forest_element_face_neighbor.cxx @@ -0,0 +1,215 @@ +/* + 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) 2024 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 + +/* We want to export the whole implementation to be callable from "C" */ +T8_EXTERN_C_BEGIN (); + +/* TODO: This function seems to be untested. + * On Nov 7 2025 i found a bug in it that would have been caught by testing + * (the face number of the element was used as the cmesh tree face number, which + * is incorrect).*/ +t8_eclass_t +t8_forest_element_neighbor_eclass (const t8_forest_t forest, const t8_locidx_t ltreeid, const t8_element_t *elem, + const int face) +{ + /* Get a pointer to the tree to read its element class */ + const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, ltreeid); + const t8_scheme *scheme = t8_forest_get_scheme (forest); + if (!scheme->element_is_root_boundary (tree_class, elem, face)) { + /* The neighbor element is inside the current tree. */ + return tree_class; + } + /* We now compute the cmesh face neighbor class across the corresponding face. + * To do so, we first need to compute the face of the root tree corresponding to the + * face of the element. */ + const int tree_face = scheme->element_get_tree_face (tree_class, elem, face); + // Debug check if tree_face is a valid face of the tree. + // Note that if elem would not be a boundary element, the return value of + // element_get_tree_face could still be within these bounds, even though it does not make sense. + // We catch that case by checking element_is_root_boundary above. + T8_ASSERT (0 <= tree_face && tree_face < t8_eclass_num_faces[tree_class]); + + const t8_locidx_t cmesh_local_tree_id = t8_forest_ltreeid_to_cmesh_ltreeid (forest, ltreeid); + const t8_cmesh_t cmesh = t8_forest_get_cmesh (forest); + return t8_cmesh_get_tree_face_neighbor_eclass (cmesh, cmesh_local_tree_id, tree_face); +} + +/* TODO: If the forest has no ghosts, then skip the ghosts + parts. In that case, process boundary elements will have 0 neighbors. +*/ +t8_gloidx_t +t8_forest_element_face_neighbor (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *elem, t8_element_t *neigh, + const t8_eclass_t neigh_eclass, int face, int *neigh_face) +{ + /* Get a pointer to the tree to read its element class */ + const t8_eclass_t eclass = t8_forest_get_tree_class (forest, ltreeid); + const t8_scheme *scheme = t8_forest_get_scheme (forest); + if (neigh_eclass == eclass && scheme->element_get_face_neighbor_inside (eclass, elem, neigh, face, neigh_face)) { + /* The neighbor was constructed and is inside the current tree. */ + return t8_forest_global_tree_id (forest, ltreeid); + } + else { + /* The neighbor does not lie inside the current tree. The content of neigh is undefined right now. */ + t8_element_t *face_element; + int tree_neigh_face; + int orientation; + int is_smaller, eclass_compare; + + const t8_cmesh_t cmesh = forest->cmesh; + /* Get the scheme associated to the element class of the boundary element. */ + /* Compute the face of elem_tree at which the face connection is. */ + const int tree_face = scheme->element_get_tree_face (eclass, elem, face); + /* compute coarse tree id */ + const t8_locidx_t lctree_id = t8_forest_ltreeid_to_cmesh_ltreeid (forest, ltreeid); + T8_ASSERT (lctree_id >= 0); +#if T8_ENABLE_DEBUG + const bool cmesh_tree_is_local = t8_cmesh_treeid_is_local_tree (cmesh, lctree_id); + T8_ASSERT (cmesh_tree_is_local || t8_cmesh_treeid_is_ghost (cmesh, lctree_id)); +#endif + + const t8_locidx_t neighbor_ctreeid + = t8_cmesh_get_face_neighbor (cmesh, lctree_id, tree_face, &tree_neigh_face, &orientation); + + if (neighbor_ctreeid < 0) { + /* This face is a domain boundary. We do not need to continue */ + return -1; + } + + /* Get the eclass for the boundary */ + const t8_eclass_t boundary_class = (t8_eclass_t) t8_eclass_face_types[eclass][tree_face]; + /* Allocate the face element */ + scheme->element_new (boundary_class, 1, &face_element); + /* Compute the face element. */ + scheme->element_get_boundary_face (eclass, elem, face, face_element); + + /* We need to find out which face is the smaller one that is the one + * according to which the orientation was computed. + * face_a is smaller then face_b if either eclass_a < eclass_b + * or eclass_a = eclass_b and face_a < face_b. */ + /* -1 eclass < neigh_eclass, 0 eclass = neigh_eclass, 1 eclass > neigh_eclass */ + eclass_compare = t8_eclass_compare (eclass, neigh_eclass); + is_smaller = 0; + if (eclass_compare == -1) { + /* The face in the current tree is the smaller one */ + is_smaller = 1; + } + else if (eclass_compare == 1) { + /* The face in the other tree is the smaller one */ + is_smaller = 0; + } + else { + + T8_ASSERT (eclass_compare == 0); + /* Check if the face of the current tree has a smaller index then the face of the neighbor tree. */ + is_smaller = tree_face <= tree_neigh_face; + } + + /* We now transform the face element to the other tree. */ + const int sign + = t8_eclass_face_orientation[eclass][tree_face] == t8_eclass_face_orientation[neigh_eclass][tree_neigh_face]; + scheme->element_transform_face (boundary_class, face_element, face_element, orientation, sign, is_smaller); + /* And now we extrude the face to the new neighbor element */ + *neigh_face = scheme->element_extrude_face (neigh_eclass, face_element, neigh, tree_neigh_face); + /* Free the face_element */ + scheme->element_destroy (boundary_class, 1, &face_element); + + const t8_gloidx_t global_neigh_id = t8_cmesh_get_global_id (cmesh, neighbor_ctreeid); + + return global_neigh_id; + } +} + +/* This function is declared in t8_forest_private.h */ +t8_gloidx_t +t8_forest_element_half_face_neighbors (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *elem, + t8_element_t *neighs[], t8_eclass_t neigh_class, int face, int num_neighs, + int dual_faces[]) +{ + if (num_neighs <= 0) { + // There are no face neighbors. + // This case might happen and we need to catch it here + // before we use neighs[0] which might not be allocated. + return -1; + } + // Use the first allocated element temporarily as same level neighbors. + t8_element_t *same_level_neighbor = neighs[0]; + + // Compute the same level neighbor element. + int same_level_dual_face; + const t8_gloidx_t neighbor_tree = t8_forest_element_face_neighbor (forest, ltreeid, elem, same_level_neighbor, + neigh_class, face, &same_level_dual_face); + + if (neighbor_tree < 0) { + // No face neighbor exists + return -1; + } + // Double check that there is a neighbor tree. + // We expect the user to only call this function if neighbors exists (neigh_calls being an input argument). + T8_ASSERT (neighbor_tree >= 0 && neighbor_tree < t8_forest_get_num_global_trees (forest)); + + /* The scheme for the current forest */ + const t8_scheme *scheme = t8_forest_get_scheme (forest); + // Check that we are allowed to refine the neighbor element. + SC_CHECK_ABORT (scheme->element_get_level (neigh_class, same_level_neighbor) < t8_forest_get_maxlevel (forest), + "Trying to refine an element beyond its maximum allowed level."); + // Double check the number of neighbors. + T8_ASSERT (num_neighs + == scheme->element_get_num_face_children (neigh_class, same_level_neighbor, same_level_dual_face)); + + // Build the half face neighbors by constructing the children at the face. + scheme->element_get_children_at_face (neigh_class, same_level_neighbor, same_level_dual_face, neighs, num_neighs, + NULL); + + // We now need to compute the dual faces of the children. + // We do this with the scheme function + if (dual_faces != NULL) { + for (int iface_child = 0; iface_child < num_neighs; ++iface_child) { + dual_faces[iface_child] + = scheme->element_face_get_child_face (neigh_class, same_level_neighbor, same_level_dual_face, iface_child); + } + } + + return neighbor_tree; +} + +int +t8_forest_leaf_face_orientation (t8_forest_t forest, const t8_locidx_t ltreeid, const t8_scheme *scheme, + const t8_element_t *leaf, int face) +{ + int orientation = 0; + const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, ltreeid); + if (scheme->element_is_root_boundary (tree_class, leaf, face)) { + t8_cmesh_t cmesh = t8_forest_get_cmesh (forest); + t8_locidx_t ltreeid_in_cmesh = t8_forest_ltreeid_to_cmesh_ltreeid (forest, ltreeid); + int iface_in_tree = scheme->element_get_tree_face (tree_class, leaf, face); + t8_cmesh_get_face_neighbor (cmesh, ltreeid_in_cmesh, iface_in_tree, NULL, &orientation); + } + + return orientation; +} + +T8_EXTERN_C_END (); diff --git a/src/t8_forest/t8_forest_neighbor/t8_forest_element_face_neighbor.h b/src/t8_forest/t8_forest_neighbor/t8_forest_element_face_neighbor.h new file mode 100644 index 0000000000..2c08609569 --- /dev/null +++ b/src/t8_forest/t8_forest_neighbor/t8_forest_element_face_neighbor.h @@ -0,0 +1,91 @@ +/* + 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) 2015 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** \file t8_forest_element_face_neighbor.h + * Functions to query for the same level face neighbors of a forest element or ghost. + * \ref t8_forest_element_neighbor_eclass + * \ref t8_forest_element_face_neighbor + * \ref t8_forest_leaf_face_orientation + */ + +#ifndef T8_FOREST_ELEMENT_FACE_NEIGHBOR_H +#define T8_FOREST_ELEMENT_FACE_NEIGHBOR_H + +#include +#include + +T8_EXTERN_C_BEGIN (); + +/** Return the eclass of the tree in which a face neighbor of a given element or ghost + * lies. + * \param [in] forest A committed forest. + * \param [in] ltreeid The local tree or ghost tree in which the element lies. 0 <= \a ltreeid < num_local_trees + num_ghost_trees + * \param [in] elem An element or ghost in the tree \a ltreeid. + * \param [in] face A face number of \a elem. + * \return The eclass of the local tree or ghost tree that + * is face neighbor of \a elem across \a face. + * T8_ECLASS_INVALID if no neighbor exists. + */ +t8_eclass_t +t8_forest_element_neighbor_eclass (const t8_forest_t forest, const t8_locidx_t ltreeid, const t8_element_t *elem, + const int face); + +/** Construct the face neighbor of an element, possibly across tree boundaries. + * Returns the global tree-id of the tree in which the neighbor element lies in. + * + * \param [in] forest The forest. + * \param [in] ltreeid The local tree in which the element lies. + * \param [in] elem The element to be considered. + * \param [in,out] neigh On input an allocated element of the scheme of the + * face_neighbors eclass. + * On output, this element's data is filled with the + * data of the face neighbor. If the neighbor does not exist + * the data could be modified arbitrarily. + * \param [in] neigh_eclass The eclass of \a neigh. + * \param [in] face The number of the face along which the neighbor should be constructed. + * \param [out] neigh_face The number of the face viewed from perspective of \a neigh. + * \return The global tree-id of the tree in which \a neigh is in. + * -1 if there exists no neighbor across that face. Domain boundary. + * -2 if the neighbor is not in a local tree or ghost tree. Process/Ghost boundary. + */ +t8_gloidx_t +t8_forest_element_face_neighbor (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *elem, t8_element_t *neigh, + const t8_eclass_t neigh_eclass, int face, int *neigh_face); + +/** Compute the leaf face orientation at given face in a forest. + * \param [in] forest The forest. Must have a valid ghost layer. + * \param [in] ltreeid A local tree id. + * \param [in] scheme The eclass scheme of the element. + * \param [in] leaf A leaf in tree \a ltreeid of \a forest. + * \param [in] face The index of the face across which the face neighbors + * are searched. + * \return Face orientation encoded as integer. + * + * For more information about the encoding of face orientation refer to \ref t8_cmesh_get_face_neighbor. + */ +int +t8_forest_leaf_face_orientation (t8_forest_t forest, const t8_locidx_t ltreeid, const t8_scheme_c *scheme, + const t8_element_t *leaf, const int face); + +T8_EXTERN_C_END (); + +#endif /* !T8_FOREST_ELEMENT_FACE_NEIGHBOR_H */ diff --git a/src/t8_forest/t8_forest_neighbor/t8_forest_leaf_face_neighbor.cxx b/src/t8_forest/t8_forest_neighbor/t8_forest_leaf_face_neighbor.cxx new file mode 100644 index 0000000000..7dcc2b301f --- /dev/null +++ b/src/t8_forest/t8_forest_neighbor/t8_forest_leaf_face_neighbor.cxx @@ -0,0 +1,558 @@ +/* + 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) 2024 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 + +/* We want to export the whole implementation to be callable from "C" */ +T8_EXTERN_C_BEGIN (); + +/** Internal data used for t8_forest_leaf_face_neighbors_iterate + * to buffer face neighbor information during leaf face neighbor search + * \ref t8_forest_leaf_face_neighbors_ext + * + * Given an element and a face, the search iterates through all leaves + * at that face and stores their information in this buffer. + * After the search the entries of the buffer are used and the + * search can start again with a clean buffer. +*/ +struct t8_lfn_user_data +{ + std::vector element_indices; /**< Element indices of the found neighbors. */ + std::vector dual_faces; /**< Dual faces of the found neighbors. */ + std::vector neighbors; /**< Pointers to the actual neighbor elements. */ +}; + +static int +t8_forest_leaf_face_neighbors_iterate (const t8_forest_t forest, const t8_locidx_t ltreeid, + [[maybe_unused]] const t8_element_t *element, const int face, const int is_leaf, + [[maybe_unused]] const t8_element_array_t *const leaf_elements, + const t8_locidx_t tree_leaf_index, void *user_data) +{ + // Output of iterate_faces: + // Array of indices in tree_leaves of all the face neighbor elements + // Assign pneighbor_leaves + // Assign dual_faces + // Assign pelement_indices + if (!is_leaf) { + // continue search until leaf level + return 1; + } + T8_ASSERT (is_leaf); + // Query whether this tree is a ghost and if so + // compute its id as a ghost tree ( 0 <= id < num_ghost_trees) + const bool is_ghost_tree = !t8_forest_tree_is_local (forest, ltreeid); + const t8_locidx_t adjusted_tree_id = !is_ghost_tree ? ltreeid : ltreeid - t8_forest_get_num_local_trees (forest); + T8_ASSERT (t8_forest_element_is_leaf_or_ghost (forest, element, adjusted_tree_id, is_ghost_tree)); + + struct t8_lfn_user_data *lfn_data = reinterpret_cast (user_data); + // face is the face of the considered leaf neighbor element and thus the + // corresponding dual face + t8_debugf ("Adding new face neighbor (leaf index %i) with dual face %i.\n", tree_leaf_index, face); + lfn_data->dual_faces.push_back (face); + // Compute the index of the element + const t8_locidx_t num_local_elements = t8_forest_get_local_num_leaf_elements (forest); + const t8_locidx_t tree_offset + = !is_ghost_tree ? t8_forest_get_tree_element_offset (forest, ltreeid) + : t8_forest_ghost_get_tree_element_offset (forest, adjusted_tree_id) + num_local_elements; + const t8_locidx_t element_index = tree_offset + tree_leaf_index; + lfn_data->element_indices.push_back (element_index); + // Add the pointer to the current element + const t8_element_t *&pnew_element = lfn_data->neighbors.emplace_back (); + if (!is_ghost_tree) { + pnew_element = t8_forest_get_leaf_element_in_tree (forest, ltreeid, tree_leaf_index); + } + else { + pnew_element = t8_forest_ghost_get_leaf_element (forest, adjusted_tree_id, tree_leaf_index); + } + return 1; +} + +/* + Set the proper return values of the leaf face neighbor computation + in case that no neighbors are found. + */ +static void +t8_forest_leaf_face_neighbors_set_no_neighbor_return_value (const t8_element_t **pneighbor_leaves[], int *dual_faces[], + int *num_neighbors, t8_locidx_t **pelement_indices, + t8_gloidx_t *gneigh_tree) +{ + *dual_faces = NULL; + *num_neighbors = 0; + *pelement_indices = NULL; + if (pneighbor_leaves) { + *pneighbor_leaves = NULL; + } + if (gneigh_tree != NULL) { + *gneigh_tree = -1; + } +} + +/* TODO: If the forest has no ghosts, then skip the ghosts + parts. In that case, process boundary elements will have 0 neighbors. +*/ +void +t8_forest_leaf_face_neighbors_ext (const t8_forest_t forest, const t8_locidx_t ltreeid, + const t8_element_t *leaf_or_ghost, const t8_element_t **pneighbor_leaves[], + const int face, int *dual_faces[], int *num_neighbors, + t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass, t8_gloidx_t *gneigh_tree, + int *orientation) +{ + /* We compute all face neighbor leaf elements of E via the following strategy: + * - Compute the same level face neighbor N + * - The neighbor tree could be a local tree or ghost (or both), + * for each variant get the leaf array of the neighbor tree and search in it: + * - Search for the first leaf element L overlapping with N. + * If it exists, it is either an ancestor or descendant of N (or N itself which is included in both definitions). + * If it does not exist, there are not leaf face neighbors in this tree. + * - If L is an ancestor of N (i.e. Level(L) <= Level(N)) then L is the only face neighbor. + * - Otherwise (Level(L) > Level (N)) we use a recursive face search across N's neighbor face, + * adding all leaf elements on the face to the face neighbors. + * This search will require L (more precise its position in the tree leaf array) as input. + **/ + + T8_ASSERT (t8_forest_is_committed (forest)); + T8_ASSERT (pelement_indices != NULL); + T8_ASSERT (dual_faces != NULL); + T8_ASSERT (num_neighbors != NULL); + +#if T8_ENABLE_DEBUG + const bool tree_is_local = t8_forest_tree_is_local (forest, ltreeid); + if (tree_is_local) { + T8_ASSERT (t8_forest_element_is_leaf (forest, leaf_or_ghost, ltreeid)); + } + else { + const t8_locidx_t local_ghost_treeid = ltreeid - t8_forest_get_num_local_trees (forest); + T8_ASSERT (t8_forest_element_is_ghost (forest, leaf_or_ghost, local_ghost_treeid)); + } +#endif + SC_CHECK_ABORT (!forest->incomplete_trees, "Leaf face neighbor is not supported for " + "forests with deleted elements.\n"); + + const t8_scheme *scheme = t8_forest_get_scheme (forest); + + if (orientation) { + // Compute the orientation of the face neighbor connection + *orientation = t8_forest_leaf_face_orientation (forest, ltreeid, scheme, leaf_or_ghost, face); + } + + /* At first we compute the same lave face neighbor element of leaf. For this, we need the + * neighbor tree's eclass and scheme. */ + const t8_eclass_t neigh_class = t8_forest_element_neighbor_eclass (forest, ltreeid, leaf_or_ghost, face); + if (neigh_class == T8_ECLASS_INVALID) { + // No face neighbor exists. + t8_forest_leaf_face_neighbors_set_no_neighbor_return_value (pneighbor_leaves, dual_faces, num_neighbors, + pelement_indices, gneigh_tree); + return; + } + if (pneigh_eclass != NULL) { + *pneigh_eclass = neigh_class; + } + + // Compute the same level face neighbor + t8_element_t *same_level_neighbor; + scheme->element_new (neigh_class, 1, &same_level_neighbor); + int same_level_neighbor_dual_face; + const t8_gloidx_t computed_gneigh_tree = t8_forest_element_face_neighbor ( + forest, ltreeid, leaf_or_ghost, same_level_neighbor, neigh_class, face, &same_level_neighbor_dual_face); + t8_debugf ("Computed same level neighbor with dual face %i\n", same_level_neighbor_dual_face); + + if (computed_gneigh_tree < 0) { + // There is no face neighbor across this face + scheme->element_destroy (neigh_class, 1, &same_level_neighbor); + t8_forest_leaf_face_neighbors_set_no_neighbor_return_value (pneighbor_leaves, dual_faces, num_neighbors, + pelement_indices, gneigh_tree); + return; + } + + // The neighbor leaves could be distributed across a local tree and a ghost + // tree. We thus possibly need to search in two different arrays. + // We store these in a vector and iterate over the entries. + // The leaf arrays themself do not store any information about their tree, + // whether it is local or ghost. + // We thus need to add this info and hence store a pair of element array and + // a bool that is true if and only if the element array corresponds to a ghost tree. + using neighbor_leaf_array = std::pair; + + // We compute the owners of the first and last face descendant. + // If the current rank is in between then the local process might have neighbor elements + // and we search the local tree. + // If other processes are in the interval of owners (lower_bound < q != p < upper_bound), + // then we (additionally or alone) search the ghost tree. + int face_owners_lower_bound = 0; + int face_owners_upper_bound = forest->mpisize - 1; + const int mpirank = forest->mpirank; + t8_forest_element_owners_at_face_bounds (forest, computed_gneigh_tree, same_level_neighbor, neigh_class, + same_level_neighbor_dual_face, &face_owners_lower_bound, + &face_owners_upper_bound); + + std::vector leaf_arrays; + + const t8_locidx_t local_neighbor_tree = t8_forest_get_local_id (forest, computed_gneigh_tree); + if (face_owners_lower_bound <= mpirank && mpirank <= face_owners_upper_bound) { + // Add the local neighbor tree's elements to the search array. + // Compute the local id of the neighbor tree and check if it is a local tree + t8_debugf ("Adding local tree to search.\n"); + if (0 <= local_neighbor_tree) { + // The neighbor tree is a local tree and hence there may be local neighbor elements. + const t8_element_array_t *tree_leaves = t8_forest_tree_get_leaf_elements (forest, local_neighbor_tree); + if (tree_leaves != nullptr) { + neighbor_leaf_array *leaf_array = new neighbor_leaf_array (tree_leaves, false); + leaf_arrays.push_back (leaf_array); + } + } + } + + if (forest->ghosts != NULL) { + if (face_owners_lower_bound != mpirank || face_owners_upper_bound != mpirank) { + // Add the neighbor tree ghost elements to the search array + t8_debugf ("Adding ghost tree to search.\n"); + const t8_locidx_t local_neighbor_ghost_treeid = t8_forest_ghost_get_ghost_treeid (forest, computed_gneigh_tree); + if (local_neighbor_ghost_treeid >= 0) { + // The neighbor tree is also a ghost tree and face neighbors of our element might + // be ghost elements. + // We add the ghost elements of that tree to our search array. + const t8_element_array_t *ghost_leaves + = t8_forest_ghost_get_tree_leaf_elements (forest, local_neighbor_ghost_treeid); + if (ghost_leaves != nullptr) { + neighbor_leaf_array *leaf_array = new neighbor_leaf_array (ghost_leaves, true); + leaf_arrays.push_back (leaf_array); + } + } + } + } + + struct t8_lfn_user_data user_data; + + // Now we iterate over the leaf arrays of the neighbor tree + // or neighbor ghost tree and find all leaf face neighbors of the element. + *num_neighbors = 0; + // Since we use REALLOC later to allocate memory of the following + // three pointers, we have to set them to NULL manually. + // This will trigger REALLOC to allocate the memory in the initial call. + // Not setting them to NULL but keeping them possibly uninitialized, will + // call REALLOC on uninitialized memory and result in memory errors. + if (pneighbor_leaves != NULL) { + /* Only set *pneighbor_leaves if a computation is desired. */ + *pneighbor_leaves = NULL; + } + + *pelement_indices = NULL; + *dual_faces = NULL; + for (auto &leaf_array : leaf_arrays) { + auto &tree_leaves = leaf_array->first; + const bool leaf_array_is_ghost = leaf_array->second; + T8_ASSERT (tree_leaves != NULL); + const t8_element_t *first_descendant; + /* + * Compute the index of the first leaf in tree_leaves that is an ancestor or descendant of + * the same_level_neighbor (might be the neighbor itself). + * Such an element might not exist in which case there are no neighbors in this tree_leaves + * array. + */ + const t8_locidx_t first_leaf_index + = t8_forest_bin_search_first_descendant_ancestor (tree_leaves, same_level_neighbor, &first_descendant); + + if (first_leaf_index >= 0) { + T8_ASSERT (first_descendant != nullptr); + const int neighbor_level = scheme->element_get_level (neigh_class, same_level_neighbor); + const int first_desc_level = scheme->element_get_level (neigh_class, first_descendant); + /* If the same level neighbor is coarser than the first found leaf, then + * we iterate over the faces of the same level neighbor. + * Otherwise, there is only one face neighbor, the first_descendant. + * We will do the iteration over the first_descendant nevertheless, but it will stop immediately. + */ + T8_ASSERT (neighbor_level >= 0); + T8_ASSERT (first_desc_level >= 0); + const bool neighbor_unique = first_desc_level <= neighbor_level; + const t8_element_t *search_this_element = neighbor_unique ? first_descendant : same_level_neighbor; + t8_debugf ("[H] Starting face search. neigh level %i, desc level %i\n", neighbor_level, first_desc_level); + + int temp_dual_face; + if (neighbor_unique && first_desc_level <= neighbor_level) { + temp_dual_face = scheme->element_face_get_ancestor_face (neigh_class, same_level_neighbor, first_desc_level, + same_level_neighbor_dual_face); + } // end if neighbor_unique + + const int search_element_dual_face = neighbor_unique ? temp_dual_face : same_level_neighbor_dual_face; + t8_debugf ("\tSearch dual face is %i\n", search_element_dual_face); + + // There may be face neighbors in this leaf array. + + // We need to restrict the array such that it contains only elements inside the search element. + // Thus, we create a new view containing all elements starting at first_leaf_index. + t8_element_array_t reduced_leaves; + if (neighbor_unique) { + t8_debugf ("Starting search with element indices %i to %i (including).\n", first_leaf_index, first_leaf_index); + t8_element_array_init_view (&reduced_leaves, tree_leaves, first_leaf_index, 1); + } + else { + /* We need to compute the first element that is not longer contained in the same_level_neighbor. + * To do so, we compute the successor of the same_level_neighbor and do + * an upper search for it in the leaf array. + * The found element (if existing) is the first leaf that is not a descendant of same_level_neighbor. */ + /* The successor might not exist because the same level neighbor is the last + * element of its level in the tree. + * To identify this case, we check whether the last leaf of the tree is an + * ancestor of same_level_neighbor. If so, then it is automatically the last leaf + * that we need to check. + * If not, we build the successor of same_level_neighbor. */ + const t8_locidx_t leaf_count = t8_element_array_get_count (tree_leaves); + const t8_element_t *last_leaf = t8_element_array_index_locidx (tree_leaves, leaf_count - 1); + T8_ASSERT (last_leaf != NULL); + t8_locidx_t last_search_element_index = -1; + if (scheme->element_is_ancestor (neigh_class, same_level_neighbor, last_leaf)) { + last_search_element_index = leaf_count - 1; + } + else { + + t8_element_t *successor; + scheme->element_new (neigh_class, 1, &successor); + scheme->element_construct_successor (neigh_class, same_level_neighbor, successor); + const int successor_level = scheme->element_get_level (neigh_class, successor); + const t8_linearidx_t successor_id = scheme->element_get_linear_id (neigh_class, successor, successor_level); + scheme->element_destroy (neigh_class, 1, &successor); + const t8_locidx_t upper_search_for_successor + = t8_forest_bin_search_upper (tree_leaves, successor_id, successor_level); + // The first index of a non descendant is the found element or the end of the array + // if no element was found. + // The last index in our search range is 1 less. + last_search_element_index = upper_search_for_successor < 0 ? leaf_count - 1 : upper_search_for_successor - 1; + } + const size_t reduced_leaf_count = last_search_element_index - first_leaf_index + 1; + T8_ASSERT (reduced_leaf_count > 0); + t8_debugf ("Starting search with element indices %i to %i (including).\n", first_leaf_index, + last_search_element_index); + t8_element_array_init_view (&reduced_leaves, tree_leaves, first_leaf_index, reduced_leaf_count); + } + // Iterate over all leaves at the face and collect them as neighbors. + const t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (forest); + // Compute the local or ghost tree id depending on whether this leaf array corresponds to a local + // tree or ghost tree. + const t8_locidx_t face_iterate_tree_id + = leaf_array_is_ghost ? t8_forest_ghost_get_ghost_treeid (forest, computed_gneigh_tree) + num_local_trees + : local_neighbor_tree; + t8_forest_iterate_faces (forest, face_iterate_tree_id, search_this_element, search_element_dual_face, + &reduced_leaves, first_leaf_index, t8_forest_leaf_face_neighbors_iterate, &user_data); + // Output of iterate_faces: + // Array of indices in tree_leaves of all the face neighbor elements + // Assign pneighbor_leaves + // Assign dual_faces + // Assign pelement_indices + // (all as growing std::vectors, resp t8_element_array) + + // num_neighbors counts the already inserted neighbors before this tree + // num_neighbors_current_tree counts the neighbors added in this tree + // total_num_neighbors temporarily counts all inserted neighbors, including this tree + const int num_neighbors_current_tree = user_data.neighbors.size (); + const int total_num_neighbors = *num_neighbors + num_neighbors_current_tree; + t8_debugf ("Found %i neighbors in tree. Adding up to %i total neighbors.\n", num_neighbors_current_tree, + total_num_neighbors); + // Copy neighbor element pointers + if (num_neighbors_current_tree > 0) { + if (pneighbor_leaves != NULL) { + *pneighbor_leaves = T8_REALLOC (*pneighbor_leaves, const t8_element_t *, total_num_neighbors); + T8_ASSERT (*pneighbor_leaves != NULL); + // Copy the pointers to pneighbor_leaves + for (t8_locidx_t ielem = 0; ielem < num_neighbors_current_tree; ++ielem) { + (*pneighbor_leaves)[ielem] = user_data.neighbors.data ()[*num_neighbors + ielem]; + } + } + // Copy element indices + *pelement_indices = T8_REALLOC (*pelement_indices, t8_locidx_t, total_num_neighbors); + T8_ASSERT (*pelement_indices != NULL); + memcpy (*pelement_indices + *num_neighbors, user_data.element_indices.data () + *num_neighbors, + num_neighbors_current_tree * sizeof (t8_locidx_t)); + // Copy dual face + *dual_faces = T8_REALLOC (*dual_faces, int, total_num_neighbors); + T8_ASSERT (*dual_faces != NULL); + memcpy (*dual_faces + *num_neighbors, user_data.dual_faces.data () + *num_neighbors, + num_neighbors_current_tree * sizeof (int)); + *num_neighbors = total_num_neighbors; + } + } // End if neighbors exist (first_leaf_index > 0) + // clean up memory allocated with new + delete leaf_array; + } + scheme->element_destroy (neigh_class, 1, &same_level_neighbor); +#if T8_ENABLE_DEBUG + // Debugging checks + if (tree_is_local && forest->ghosts != NULL) { + // For local elements we must have found face neighbors by now. + T8_ASSERT (*num_neighbors > 0); + } + // All neighbor elements must be valid + if (pneighbor_leaves != NULL) { + for (int ineigh = 0; ineigh < *num_neighbors; ++ineigh) { + T8_ASSERT (scheme->element_is_valid (neigh_class, (*pneighbor_leaves)[ineigh])); + t8_debugf ("Face neighbor %p is valid.\n", (void *) (*pneighbor_leaves)[ineigh]); + scheme->element_debug_print (neigh_class, (*pneighbor_leaves)[ineigh]); + } + } +#endif // T8_ENABLE_DEBUG + // If no neighbors are found, set the proper return values. + if (*num_neighbors == 0) { + t8_debugf ("Found no neighbors\n"); + t8_forest_leaf_face_neighbors_set_no_neighbor_return_value (pneighbor_leaves, dual_faces, num_neighbors, + pelement_indices, gneigh_tree); + T8_ASSERT (*dual_faces == NULL); + T8_ASSERT (*pelement_indices == NULL); + T8_ASSERT (pneighbor_leaves == NULL || *pneighbor_leaves == NULL); + T8_ASSERT (gneigh_tree == NULL || *gneigh_tree == -1); + } + + if (gneigh_tree != NULL) { + *gneigh_tree = computed_gneigh_tree; + } +} + +void +t8_forest_leaf_face_neighbors (const t8_forest_t forest, const t8_locidx_t ltreeid, const t8_element_t *leaf, + const t8_element_t **pneighbor_leaves[], const int face, int *dual_faces[], + int *num_neighbors, t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass) +{ + t8_forest_leaf_face_neighbors_ext (forest, ltreeid, leaf, pneighbor_leaves, face, dual_faces, num_neighbors, + pelement_indices, pneigh_eclass, NULL, NULL); +} + +t8_locidx_t +t8_forest_same_level_leaf_face_neighbor_index (const t8_forest_t forest, const t8_locidx_t element_index, + const int face_index, const t8_gloidx_t global_treeid, int *dual_face) +{ + const t8_locidx_t num_local_elements = t8_forest_get_local_num_leaf_elements (forest); +#if T8_ENABLE_DEBUG + const t8_locidx_t num_ghosts = t8_forest_get_num_ghosts (forest); + T8_ASSERT (0 <= element_index && element_index < num_local_elements + num_ghosts); +#endif + const bool is_local = element_index < num_local_elements; + + t8_locidx_t local_tree; + t8_locidx_t element_index_in_tree; + const t8_element_t *element; + if (is_local) { + local_tree = t8_forest_get_local_id (forest, global_treeid); + element_index_in_tree = element_index - t8_forest_get_tree_element_offset (forest, local_tree); + element = t8_forest_get_leaf_element_in_tree (forest, local_tree, element_index_in_tree); + } + else { + local_tree = t8_forest_ghost_get_ghost_treeid (forest, global_treeid); + const t8_locidx_t ghost_offset_in_tree = t8_forest_ghost_get_tree_element_offset (forest, local_tree); + element_index_in_tree = element_index - num_local_elements - ghost_offset_in_tree; + element = t8_forest_ghost_get_leaf_element (forest, local_tree, element_index_in_tree); + local_tree += t8_forest_get_num_local_trees (forest); + } + + int *dual_faces; + int num_neighbors = 0; + t8_locidx_t *element_indices; + t8_eclass_t neigh_class; + + t8_debugf ("Same level leaf neighbor for index %i. Which is %s element %i in tree %i.\n", element_index, + element_index < num_local_elements ? "local" : "ghost", element_index_in_tree, local_tree); + + t8_forest_leaf_face_neighbors (forest, local_tree, element, NULL, face_index, &dual_faces, &num_neighbors, + &element_indices, &neigh_class); + + T8_ASSERT (num_neighbors == 0 || num_neighbors == 1); + + if (num_neighbors == 0) { + *dual_face = -1; + return -1; + } + + *dual_face = dual_faces[0]; + const t8_locidx_t neigh_index = element_indices[0]; + + T8_FREE (element_indices); + T8_FREE (dual_faces); + + return neigh_index; +} + +/* This function is declared in t8_forest_private.h */ +void +t8_forest_print_all_leaf_neighbors (t8_forest_t forest) +{ + t8_locidx_t ltree, ielem; + const t8_element_t **neighbor_leaves; + int iface, num_neighbors, ineigh; + t8_eclass_t eclass, neigh_eclass; + const t8_scheme *scheme = t8_forest_get_scheme (forest); + t8_locidx_t *element_indices; + int *dual_faces; + char buffer[BUFSIZ]; + int allocate_first_desc = 0, allocate_tree_offset = 0; + int allocate_el_offset = 0; + + if (forest->tree_offsets == NULL) { + allocate_tree_offset = 1; + t8_forest_partition_create_tree_offsets (forest); + } + if (forest->global_first_desc == NULL) { + allocate_first_desc = 1; + t8_forest_partition_create_first_desc (forest); + } + if (forest->element_offsets == NULL) { + allocate_el_offset = 1; + t8_forest_partition_create_offsets (forest); + } + for (ielem = 0; ielem < t8_forest_get_local_num_leaf_elements (forest); ielem++) { + /* Get a pointer to the ielem-th element, its eclass, treeid and scheme */ + const t8_element_t *leaf = t8_forest_get_leaf_element (forest, ielem, <ree); + eclass = t8_forest_get_tree_class (forest, ltree); + /* Iterate over all faces */ + for (iface = 0; iface < scheme->element_get_num_faces (eclass, leaf); iface++) { + t8_forest_leaf_face_neighbors (forest, ltree, leaf, &neighbor_leaves, iface, &dual_faces, &num_neighbors, + &element_indices, &neigh_eclass); + t8_debugf ("Element %li across face %i has %i leaf neighbors (with dual faces).\n", (long) ielem, iface, + num_neighbors); + snprintf (buffer, BUFSIZ, "\tIndices:\t"); + for (ineigh = 0; ineigh < num_neighbors; ineigh++) { + snprintf (buffer + strlen (buffer), BUFSIZ - strlen (buffer), "%li (%i) ", (long) element_indices[ineigh], + dual_faces[iface]); + } + t8_debugf ("%s\n", buffer); + if (num_neighbors > 0) { + T8_FREE (element_indices); + T8_FREE (neighbor_leaves); + T8_FREE (dual_faces); + } + } + } + if (allocate_tree_offset) { + t8_shmem_array_destroy (&forest->tree_offsets); + } + if (allocate_first_desc) { + t8_shmem_array_destroy (&forest->global_first_desc); + } + if (allocate_el_offset) { + t8_shmem_array_destroy (&forest->element_offsets); + } +} + +T8_EXTERN_C_END (); diff --git a/src/t8_forest/t8_forest_neighbor/t8_forest_leaf_face_neighbor.h b/src/t8_forest/t8_forest_neighbor/t8_forest_leaf_face_neighbor.h new file mode 100644 index 0000000000..d151936e28 --- /dev/null +++ b/src/t8_forest/t8_forest_neighbor/t8_forest_leaf_face_neighbor.h @@ -0,0 +1,131 @@ +/* + 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) 2015 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** \file t8_forest_leaf_face_neighbor.h + * Functions to query for the leaf face neighbors of a forest leaf element or leaf ghost. + * \ref t8_forest_leaf_face_neighbors + * \ref t8_forest_leaf_face_neighbors_ext + * \ref t8_forest_same_level_leaf_face_neighbor_index + */ + +#ifndef T8_FOREST_LEAF_FACE_NEIGHBOR_H +#define T8_FOREST_LEAF_FACE_NEIGHBOR_H + +#include +#include + +T8_EXTERN_C_BEGIN (); + +/** Compute the leaf face neighbors of a forest leaf element or ghost leaf. + * \param [in] forest The forest. + * \param [in] ltreeid A local tree id (could also be a ghost tree). 0 <= \a ltreeid < num_local trees+num_ghost_trees + * \param [in] leaf A leaf in tree \a ltreeid of \a forest. + * \param [out] pneighbor_leaves Unallocated on input. On output the neighbor + * leaves are stored here. + * \param [in] face The index of the face across which the face neighbors + * are searched. + * \param [out] dual_faces On output the face id's of the neighboring elements' faces. + * \param [out] num_neighbors On output the number of neighbor leaves. + * \param [out] pelement_indices Unallocated on input. On output the element indices + * of the neighbor leaves are stored here. + * 0, 1, ... num_local_el - 1 for local leaves and + * num_local_el , ... , num_local_el + num_ghosts - 1 for ghosts. + * \param [out] pneigh_eclass On output the eclass of the neighbor elements. + * \note If there are no face neighbors, then *pneighbor_leaves = NULL, num_neighbors = 0, + * and *pelement_indices = NULL on output. + * \note \a forest must be committed before calling this function. + * \note If \a forest does not have a ghost layer then leaf elements at the process boundaries have 0 neighbors. (The function output for leaf elements then depends on the parallel partition.) + * \note Important! This routine allocates memory which must be freed. Do it like this: + * + * if (num_neighbors > 0) { + * T8_FREE (pneighbor_leaves); + * T8_FREE (pelement_indices); + * T8_FREE (dual_faces); + * } + * + */ +void +t8_forest_leaf_face_neighbors (const t8_forest_t forest, const t8_locidx_t ltreeid, const t8_element_t *leaf, + const t8_element_t **pneighbor_leaves[], const int face, int *dual_faces[], + int *num_neighbors, t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass); + +/** Like \ref t8_forest_leaf_face_neighbors but also provides information about the global neighbors and the orientation. + * \param [in] forest The forest. Must have a valid ghost layer. + * \param [in] ltreeid A local tree id (could also be a ghost tree). 0 <= \a ltreeid < num_local trees+num_ghost_trees + * \param [in] leaf A leaf in tree \a ltreeid of \a forest. + * \param [out] pneighbor_leaves Unallocated on input. On output the neighbor + * leaves are stored here. + * \param [in] face The index of the face across which the face neighbors + * are searched. + * \param [out] dual_faces On output the face id's of the neighboring elements' faces. + * \param [out] num_neighbors On output the number of neighbor leaves. + * \param [out] pelement_indices Unallocated on input. On output the element indices + * of the neighbor leaves are stored here. + * 0, 1, ... num_local_el - 1 for local leaves and + * num_local_el , ... , num_local_el + num_ghosts - 1 for ghosts. + * \param [out] pneigh_eclass On output the eclass of the neighbor elements. + * \param [out] gneigh_tree The global tree IDs of the neighbor trees. + * \param [out] orientation If not NULL on input, the face orientation is computed and stored here. + * Thus, if the face connection is an inter-tree connection the orientation of the tree-to-tree connection is stored. + * Otherwise, the value 0 is stored. + * All other parameters and behavior are identical to \ref t8_forest_leaf_face_neighbors. + * \note If there are no face neighbors, then *pneighbor_leaves = NULL, num_neighbors = 0, + * and *pelement_indices = NULL on output. + * \note \a forest must be committed before calling this function. + * + * \note Important! This routine allocates memory which must be freed. Do it like this: + * + * if (num_neighbors > 0) { + * T8_FREE (pneighbor_leaves); + * T8_FREE (pelement_indices); + * T8_FREE (dual_faces); + * } + * + */ +void +t8_forest_leaf_face_neighbors_ext (const t8_forest_t forest, const t8_locidx_t ltreeid, + const t8_element_t *leaf_or_ghost, const t8_element_t **pneighbor_leaves[], + const int face, int *dual_faces[], int *num_neighbors, + t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass, t8_gloidx_t *gneigh_tree, + int *orientation); + +/** Given a leaf element or ghost index in "all local elements + ghosts" enumeration + * compute the index of the face neighbor of the element - provided that only one or no + * face neighbors exists. + * HANDLE WITH CARE. DO NOT CALL IF THE FOREST IS ADAPTED. + * + * \param[in] forest The forest. Must be committed. + * \param[in] element_index Index of an element in \a forest. Must have only one or no facen neighbors across the given face. + * 0 <= \a element_index < num_local_elements + num_ghosts + * \param[in] face_index Index of a face of \a element. + * \param[in] global_treeid Global index of the tree that contains \a element. + * \param[out] dual_face Return value, the dual_face index of the face neighbor. + * \return The index of the face neighbor leaf (local element or ghost). + * \note Do not call if you are unsure about the number of face neighbors. In particular if the forest is adapted and not uniform. + */ +t8_locidx_t +t8_forest_same_level_leaf_face_neighbor_index (const t8_forest_t forest, const t8_locidx_t element_index, + const int face_index, const t8_gloidx_t global_treeid, int *dual_face); + +T8_EXTERN_C_END (); + +#endif /* !T8_FOREST_LEAF_FACE_NEIGHBOR_H */ From fc6d695998ef01eb297094d8573c81bd19b60160 Mon Sep 17 00:00:00 2001 From: Johannes Holke Date: Tue, 18 Nov 2025 16:59:15 +0100 Subject: [PATCH 8/8] Change includes wherever face neighbors are needed --- example/advect/t8_advection.cxx | 5 +---- src/CMakeLists.txt | 4 +++- src/t8_forest/t8_forest.h | 2 ++ src/t8_forest/t8_forest_balance.cxx | 1 + src/t8_forest/t8_forest_ghost.cxx | 1 + test/t8_forest/t8_gtest_forest_face_normal.cxx | 3 +-- test/t8_forest/t8_gtest_half_neighbors.cxx | 6 +++--- tutorials/general/t8_step6_stencil.cxx | 4 +--- 8 files changed, 13 insertions(+), 13 deletions(-) diff --git a/example/advect/t8_advection.cxx b/example/advect/t8_advection.cxx index 1c7a2658d4..1f420dc45b 100644 --- a/example/advect/t8_advection.cxx +++ b/example/advect/t8_advection.cxx @@ -23,10 +23,7 @@ #include #include #include -#include -#include -#include -#include +#include #include #include #include diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 26723b22fb..fc2cb98558 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -177,7 +177,9 @@ target_sources( T8 PRIVATE t8_forest/t8_forest_ghost.cxx t8_forest/t8_forest_iterate.cxx t8_forest/t8_forest_balance.cxx - t8_forest/t8_forest_search/t8_forest_search.cxx + t8_forest/t8_forest_search/t8_forest_search.cxx + t8_forest/t8_forest_neighbor/t8_forest_element_face_neighbor.cxx + t8_forest/t8_forest_neighbor/t8_forest_leaf_face_neighbor.cxx t8_geometry/t8_geometry.cxx t8_geometry/t8_geometry_helpers.c t8_geometry/t8_geometry_base.cxx diff --git a/src/t8_forest/t8_forest.h b/src/t8_forest/t8_forest.h index 12b0f0a0c8..5f31ffcfbe 100644 --- a/src/t8_forest/t8_forest.h +++ b/src/t8_forest/t8_forest.h @@ -31,5 +31,7 @@ #include #include #include +#include +#include #endif /* !T8_FOREST_H */ diff --git a/src/t8_forest/t8_forest_balance.cxx b/src/t8_forest/t8_forest_balance.cxx index c7639c9622..48253ad0e2 100644 --- a/src/t8_forest/t8_forest_balance.cxx +++ b/src/t8_forest/t8_forest_balance.cxx @@ -27,6 +27,7 @@ #include #include #include +#include #include /* We want to export the whole implementation to be callable from "C" */ diff --git a/src/t8_forest/t8_forest_ghost.cxx b/src/t8_forest/t8_forest_ghost.cxx index e29e43120e..acaee67b61 100644 --- a/src/t8_forest/t8_forest_ghost.cxx +++ b/src/t8_forest/t8_forest_ghost.cxx @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include diff --git a/test/t8_forest/t8_gtest_forest_face_normal.cxx b/test/t8_forest/t8_gtest_forest_face_normal.cxx index fd83756b96..791a8f63de 100644 --- a/test/t8_forest/t8_gtest_forest_face_normal.cxx +++ b/test/t8_forest/t8_gtest_forest_face_normal.cxx @@ -25,8 +25,7 @@ #include #include #include -#include -#include +#include #include #include diff --git a/test/t8_forest/t8_gtest_half_neighbors.cxx b/test/t8_forest/t8_gtest_half_neighbors.cxx index 494022a48b..05b6820253 100644 --- a/test/t8_forest/t8_gtest_half_neighbors.cxx +++ b/test/t8_forest/t8_gtest_half_neighbors.cxx @@ -26,13 +26,13 @@ #include #include -#include #include +#include +#include +#include #include #include #include -#include -#include #include class forest_half_neighbors: public testing::TestWithParam, int>> { diff --git a/tutorials/general/t8_step6_stencil.cxx b/tutorials/general/t8_step6_stencil.cxx index 13cf80f3c1..99ad4c9d7d 100644 --- a/tutorials/general/t8_step6_stencil.cxx +++ b/tutorials/general/t8_step6_stencil.cxx @@ -42,9 +42,7 @@ #include /* General t8code header, always include this. */ #include /* Basic operations on 3D vectors. */ #include /* cmesh definition and basic interface. */ -#include /* forest definition and basic interface. */ -#include /* save forest */ -#include /* geometrical information of the forest */ +#include /* forest definition and basic interface. */ #include /* A collection of exemplary cmeshes */ #include /* default refinement scheme. */ #include