diff --git a/NEWS.md b/NEWS.md index 017db53cd0..0dab757cb9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -21,6 +21,24 @@ No! If your code is only a couple of lines long AND has very little impact on th # User Updates for the upcoming t8code release (February 2026 - version format unclear) +## Updates to t8_forest_leaf_face_neighbors + +The t8_forest_leaf_face_neighbors and t8_forest_leaf_face_neighbors_ext functions (see t8_forest_general.h) have undergone a major update. + +- The forest does not need to be balanced anymore to compute the leaf face neighbors. Thus, there can be any arbitrary number of leaf face neighbors for a given element. +The last argument (forest_is_balanced) was removed. + +- The allocation behaviour changed. element_destroy should not be called anymore. +The function now returns pointers to t8code internal elements in the neighbor_leaves array. +To free the memory use +```C++ +if (num_neighbors > 0) { + T8_FREE (pneighbor_leaves); + T8_FREE (pelement_indices); + T8_FREE (dual_faces); +} +``` + ## Using structs instead of classes We decided to use `struct` instead of `class` throughout the entire codebase. diff --git a/example/advect/t8_advection.cxx b/example/advect/t8_advection.cxx index 3fcdd4a07d..49a5a1ed1c 100644 --- a/example/advect/t8_advection.cxx +++ b/example/advect/t8_advection.cxx @@ -984,7 +984,7 @@ t8_advect_problem_init_elements (t8_advect_problem_t *problem) { t8_locidx_t itree, ielement, idata; t8_locidx_t num_trees, num_elems_in_tree; - t8_element_t **neighbors; + const t8_element_t **neighbors; int iface, ineigh; t8_advect_element_data_t *elem_data; const t8_scheme *scheme = t8_forest_get_scheme (problem->forest); @@ -1035,13 +1035,12 @@ t8_advect_problem_init_elements (t8_advect_problem_t *problem) t8_forest_leaf_face_neighbors (problem->forest, itree, element, &neighbors, iface, &elem_data->dual_faces[iface], &elem_data->num_neighbors[iface], - &elem_data->neighs[iface], &neigh_eclass, 1); + &elem_data->neighs[iface], &neigh_eclass); for (ineigh = 0; ineigh < elem_data->num_neighbors[iface]; ineigh++) { elem_data->neigh_level[iface] = scheme->element_get_level (neigh_eclass, neighbors[ineigh]); } if (elem_data->num_neighbors[iface] > 0) { - scheme->element_destroy (neigh_eclass, elem_data->num_neighbors[iface], neighbors); T8_FREE (neighbors); //t8_global_essentialf("alloc face %i of elem %i\n", iface, ielement); elem_data->fluxes[iface] = T8_ALLOC (double, elem_data->num_neighbors[iface]); @@ -1197,7 +1196,7 @@ t8_advect_solve (t8_cmesh_t cmesh, t8_flow_function_3d_fn u, t8_example_level_se int done = 0; int adapted_or_partitioned = 0; int dual_face; - t8_element_t **neighs; + const t8_element_t **neighs; t8_eclass_t neigh_eclass; double total_time, solve_time = 0; double ghost_exchange_time, ghost_waittime, neighbor_time, flux_time; @@ -1291,15 +1290,13 @@ t8_advect_solve (t8_cmesh_t cmesh, t8_flow_function_3d_fn u, t8_example_level_se neighbor_time = -sc_MPI_Wtime (); t8_forest_leaf_face_neighbors (problem->forest, itree, elem, &neighs, iface, &elem_data->dual_faces[iface], &elem_data->num_neighbors[iface], - &elem_data->neighs[iface], &neigh_eclass, 1); + &elem_data->neighs[iface], &neigh_eclass); for (ineigh = 0; ineigh < elem_data->num_neighbors[iface]; ineigh++) { elem_data->neigh_level[iface] = scheme->element_get_level (neigh_eclass, neighs[ineigh]); } T8_ASSERT (neighs != NULL || elem_data->num_neighbors[iface] == 0); if (neighs != NULL) { - scheme->element_destroy (neigh_eclass, elem_data->num_neighbors[iface], neighs); - T8_FREE (neighs); } diff --git a/example/forest/t8_test_face_iterate.cxx b/example/forest/t8_test_face_iterate.cxx index fe51e449c3..a1a360387c 100644 --- a/example/forest/t8_test_face_iterate.cxx +++ b/example/forest/t8_test_face_iterate.cxx @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -41,17 +42,17 @@ typedef struct } t8_test_fiterate_udata_t; static int -t8_test_fiterate_callback (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, int face, - void *user_data, t8_locidx_t leaf_index) +t8_test_fiterate_callback (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *element, int face, int is_leaf, + [[maybe_unused]] const t8_element_array_t *leaf_elements, t8_locidx_t leaf_index, + void *user_data) { - double *coords; - - if (leaf_index >= 0) { - coords = ((t8_test_fiterate_udata_t *) user_data)->coords; + if (is_leaf) { + t8_test_fiterate_udata_t *test_user_data = (t8_test_fiterate_udata_t *) user_data; + double *coords = test_user_data->coords; t8_forest_element_coordinate (forest, ltreeid, element, 0, coords); t8_debugf ("Leaf element in tree %i at face %i, tree local index %i has corner 0 coords %lf %lf %lf\n", ltreeid, face, (int) leaf_index, coords[0], coords[1], coords[2]); - ((t8_test_fiterate_udata_t *) user_data)->count++; + test_user_data->count++; } return 1; } @@ -75,27 +76,37 @@ t8_basic_adapt ([[maybe_unused]] t8_forest_t forest, [[maybe_unused]] t8_forest_ static void t8_test_fiterate (t8_forest_t forest) { - t8_locidx_t itree, num_trees; - t8_eclass_t eclass; - const t8_scheme *scheme = t8_forest_get_scheme (forest); - t8_element_t *nca; - t8_element_array_t *leaf_elements; t8_test_fiterate_udata_t udata; - int iface; - num_trees = t8_forest_get_num_local_trees (forest); - for (itree = 0; itree < num_trees; itree++) { - eclass = t8_forest_get_tree_class (forest, itree); - const t8_element_t *first_el = t8_forest_get_leaf_element_in_tree (forest, itree, 0); + t8_forest_set_user_data (forest, &udata); + const t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (forest); + const t8_locidx_t num_trees = num_local_trees + t8_forest_get_num_ghost_trees (forest); + for (t8_locidx_t itree = 0; itree < num_trees; itree++) { + const t8_eclass_t eclass = t8_forest_get_tree_class (forest, itree); + const t8_scheme *scheme = t8_forest_get_scheme (forest); + // Query whether this tree is a ghost and compute its ghost tree id. + const bool is_ghost = itree >= num_local_trees; + const t8_locidx_t ghost_tree_id = is_ghost ? itree - num_local_trees : -1; + // Get the number of elements or ghosts in this tree + const t8_locidx_t num_tree_elements = is_ghost ? t8_forest_ghost_tree_num_leaf_elements (forest, ghost_tree_id) + : t8_forest_get_tree_num_leaf_elements (forest, itree); + // Get all leaf elements + const t8_element_array_t *leaf_elements = is_ghost + ? t8_forest_ghost_get_tree_leaf_elements (forest, ghost_tree_id) + : t8_forest_tree_get_leaf_elements (forest, itree); + // Get the first and last element + const t8_element_t *first_el = (const t8_element_t *) t8_element_array_index_locidx (leaf_elements, 0); const t8_element_t *last_el - = t8_forest_get_leaf_element_in_tree (forest, itree, t8_forest_get_tree_num_leaf_elements (forest, itree) - 1); + = (const t8_element_t *) t8_element_array_index_locidx (leaf_elements, num_tree_elements - 1); + + t8_element_t *nca; scheme->element_new (eclass, 1, &nca); scheme->element_get_nca (eclass, first_el, last_el, nca); - leaf_elements = t8_forest_tree_get_leaf_elements (forest, itree); - for (iface = 0; iface < scheme->element_get_num_faces (eclass, nca); iface++) { + // + for (int iface = 0; iface < scheme->element_get_num_faces (eclass, nca); iface++) { udata.count = 0; - t8_forest_iterate_faces (forest, itree, nca, iface, leaf_elements, &udata, 0, t8_test_fiterate_callback); + t8_forest_iterate_faces (forest, itree, nca, iface, leaf_elements, 0, t8_test_fiterate_callback, &udata); t8_debugf ("Leaf elements at face %i:\t%i\n", iface, udata.count); } scheme->element_destroy (eclass, 1, &nca); @@ -139,6 +150,7 @@ t8_test_fiterate_refine_and_partition (t8_cmesh_t cmesh, int level, sc_MPI_Comm /* partition the adapted forest */ t8_forest_init (&forest_partition); t8_forest_set_partition (forest_partition, forest_adapt, 0); + t8_forest_set_ghost (forest_partition, 1, T8_GHOST_FACES); t8_forest_commit (forest_partition); t8_debugf ("Created ghost structure with %li ghost elements.\n", (long) t8_forest_get_num_ghosts (forest_partition)); if (!no_vtk) { diff --git a/mesh_handle/element.hxx b/mesh_handle/element.hxx index 78cb844f22..66c92575da 100644 --- a/mesh_handle/element.hxx +++ b/mesh_handle/element.hxx @@ -345,15 +345,14 @@ class element: public TCompetences>... { return this->m_neighbors[face]; } } - - t8_element_t** neighbors; /**< Neighboring elements. */ - int* dual_faces_internal; /**< Face indices of the neighbor elements. */ - int num_neighbors; /**< Number of neighboring elements. */ - t8_locidx_t* neighids; /**< Neighboring elements ids. */ - t8_eclass_t neigh_class; /**< Neighboring elements tree class. */ + const t8_element_t** neighbors; /**< Neighboring elements. */ + int* dual_faces_internal; /**< Face indices of the neighbor elements. */ + int num_neighbors; /**< Number of neighboring elements. */ + t8_locidx_t* neighids; /**< Neighboring elements ids. */ + t8_eclass_t neigh_class; /**< Neighboring elements tree class. */ t8_forest_leaf_face_neighbors (m_mesh->m_forest, m_tree_id, m_element, &neighbors, face, &dual_faces_internal, - &num_neighbors, &neighids, &neigh_class, t8_forest_is_balanced (m_mesh->m_forest)); + &num_neighbors, &neighids, &neigh_class); if (dual_faces) { dual_faces->get ().assign (dual_faces_internal, dual_faces_internal + num_neighbors); } @@ -369,7 +368,6 @@ class element: public TCompetences>... { } if (num_neighbors > 0) { // Free allocated memory. - t8_forest_get_scheme (m_mesh->m_forest)->element_destroy (get_tree_class (), num_neighbors, neighbors); T8_FREE (neighbors); T8_FREE (dual_faces_internal); T8_FREE (neighids); diff --git a/src/t8_forest/t8_forest.cxx b/src/t8_forest/t8_forest.cxx index 3e6993b74f..de3a86a983 100644 --- a/src/t8_forest/t8_forest.cxx +++ b/src/t8_forest/t8_forest.cxx @@ -34,6 +34,8 @@ #include #include #include +#include +#include #include #include #include @@ -49,6 +51,7 @@ #include #include +#include /* We want to export the whole implementation to be callable from "C" */ T8_EXTERN_C_BEGIN (); @@ -1422,119 +1425,84 @@ 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). + * See https://github.com/DLR-AMR/t8code/issues/2240 */ t8_eclass_t -t8_forest_element_neighbor_eclass (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *elem, int face) +t8_forest_element_neighbor_eclass (const t8_forest_t forest, const t8_locidx_t ltreeid, const t8_element_t *elem, + const int face) { - t8_ctree_t coarse_tree; - int tree_face; - t8_locidx_t lcoarse_neighbor; - t8_cmesh_t cmesh; - /* Get a pointer to the tree to read its element class */ - const t8_tree_t tree = t8_forest_get_tree (forest, ltreeid); - const t8_eclass_t tree_class = tree->eclass; + 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; } - else { - /* The neighbor is in a neighbor tree */ - /* If the face neighbor is not inside the tree, we have to find out the tree - * face and the tree's face neighbor along that face. */ - tree_face = scheme->element_get_tree_face (tree_class, elem, face); - - cmesh = t8_forest_get_cmesh (forest); - /* Get the coarse tree corresponding to tree */ - coarse_tree = t8_forest_get_coarse_tree (forest, ltreeid); - /* Get the (coarse) local id of the tree neighbor */ - lcoarse_neighbor = t8_cmesh_trees_get_face_neighbor (coarse_tree, tree_face); - T8_ASSERT (0 <= lcoarse_neighbor); - if (lcoarse_neighbor < t8_cmesh_get_num_local_trees (cmesh)) { - /* The tree neighbor is a local tree */ - return t8_cmesh_get_tree_class (cmesh, lcoarse_neighbor); - } - else { - T8_ASSERT (lcoarse_neighbor - t8_cmesh_get_num_local_trees (cmesh) < cmesh->num_ghosts); - /* The tree neighbor is a ghost */ - return t8_cmesh_get_ghost_class (cmesh, lcoarse_neighbor - t8_cmesh_get_num_local_trees (cmesh)); - } - } + /* 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: Function declaration return statement does not match the implementation. +// Check this. 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, - t8_eclass_t neigh_eclass, int face, int *neigh_face) + const t8_eclass_t neigh_eclass, int face, int *neigh_face) { /* Get a pointer to the tree to read its element class */ - const t8_tree_t tree = t8_forest_get_tree (forest, ltreeid); - const t8_eclass_t eclass = tree->eclass; + 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 ltreeid + t8_forest_get_first_local_tree_id (forest); + 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_eclass_t neigh_eclass; t8_element_t *face_element; - t8_cmesh_t cmesh; - t8_locidx_t lctree_id, lcneigh_id; - t8_locidx_t *face_neighbor; - t8_gloidx_t global_neigh_id; - t8_cghost_t ghost; - int8_t *ttf; - int tree_face, tree_neigh_face; + int tree_neigh_face; + int orientation; int is_smaller, eclass_compare; - int F, sign; - cmesh = forest->cmesh; + 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. */ - tree_face = scheme->element_get_tree_face (eclass, elem, face); + const int tree_face = scheme->element_get_tree_face (eclass, elem, face); /* compute coarse tree id */ - lctree_id = t8_forest_ltreeid_to_cmesh_ltreeid (forest, ltreeid); - if (t8_cmesh_tree_face_is_boundary (cmesh, lctree_id, tree_face)) { + const t8_locidx_t cmesh_ltreeid = 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, cmesh_ltreeid); + T8_ASSERT (cmesh_tree_is_local || t8_cmesh_treeid_is_ghost (cmesh, cmesh_ltreeid)); +#endif + + const t8_locidx_t neighbor_ctreeid + = t8_cmesh_get_face_neighbor (cmesh, cmesh_ltreeid, 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); - /* Get the coarse tree that contains elem. - * Also get the face neighbor information of the coarse tree. */ - (void) t8_cmesh_trees_get_tree_ext (cmesh->trees, lctree_id, &face_neighbor, &ttf); - /* Compute the local id of the face neighbor tree. */ - lcneigh_id = face_neighbor[tree_face]; - /* F is needed to compute the neighbor face number and the orientation. - * tree_neigh_face = ttf % F - * or = ttf / F - */ - F = t8_eclass_max_num_faces[cmesh->dimension]; - /* compute the neighbor face */ - tree_neigh_face = ttf[tree_face] % F; - if (lcneigh_id == lctree_id && tree_face == tree_neigh_face) { - /* This face is a domain boundary and there is no neighbor */ - return -1; - } - /* We now compute the eclass of the neighbor tree. */ - if (lcneigh_id < t8_cmesh_get_num_local_trees (cmesh)) { - /* The face neighbor is a local tree */ - /* Get the eclass of the neighbor tree */ - neigh_eclass = t8_cmesh_get_tree_class (cmesh, lcneigh_id); - global_neigh_id = lcneigh_id + t8_cmesh_get_first_treeid (cmesh); - } - else { - /* The face neighbor is a ghost tree */ - T8_ASSERT (cmesh->num_local_trees <= lcneigh_id && lcneigh_id < cmesh->num_ghosts + cmesh->num_local_trees); - /* Get the eclass of the neighbor tree */ - ghost = t8_cmesh_trees_get_ghost (cmesh->trees, lcneigh_id - t8_cmesh_get_num_local_trees (cmesh)); - neigh_eclass = ghost->eclass; - global_neigh_id = ghost->treeid; - } + /* 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 @@ -1556,14 +1524,18 @@ t8_forest_element_face_neighbor (t8_forest_t forest, t8_locidx_t ltreeid, const /* 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. */ - 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, ttf[tree_face] / F, sign, is_smaller); + 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; } } @@ -1636,260 +1608,476 @@ t8_forest_leaf_face_orientation (t8_forest_t forest, const t8_locidx_t ltreeid, return orientation; } -void -t8_forest_leaf_face_neighbors_ext (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *leaf, - t8_element_t **pneighbor_leaves[], int face, int *dual_faces[], int *num_neighbors, - t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass, int forest_is_balanced, - t8_gloidx_t *gneigh_tree, int *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 { - t8_eclass_t eclass; - t8_gloidx_t gneigh_treeid; - t8_locidx_t lneigh_treeid = -1; - t8_locidx_t lghost_treeid = -1, *element_indices, element_index; - const t8_scheme *scheme = t8_forest_get_scheme (forest); - const t8_element_t *ancestor; - t8_element_t **neighbor_leaves; - t8_linearidx_t neigh_id; - int num_children_at_face, at_maxlevel; - int ineigh, *owners, different_owners, have_ghosts; + 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 - t8_forest_get_num_local_trees (forest) : ltreeid; + 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_ghost_get_tree_element_offset (forest, adjusted_tree_id) + num_local_elements + : t8_forest_get_tree_element_offset (forest, ltreeid); + 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_ghost_get_leaf_element (forest, adjusted_tree_id, tree_leaf_index); + } + else { + pnew_element = t8_forest_get_leaf_element_in_tree (forest, ltreeid, tree_leaf_index); + } + return 1; +} + +/** + * Set the proper return values of the leaf face neighbor computation + * in case that no neighbors are found. + * + * \param [out] pneighbor_leaves If not NULL, will be set to NULL. + * \param [out] dual_faces Will be set to NULL. + * \param [out] num_neighbors Will be set to 0. + * \param [out] pelement_indices Will be set to NULL. + * \param [out] gneigh_tree If not NULL, will be set to -1. + */ +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; + } +} + +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 (t8_forest_element_is_leaf (forest, leaf, ltreeid)); - T8_ASSERT (!forest_is_balanced || t8_forest_is_balanced (forest)); - SC_CHECK_ABORT (forest_is_balanced, "leaf face neighbors is not implemented " - "for unbalanced forests.\n"); /* TODO: write version for unbalanced forests */ - SC_CHECK_ABORT (forest->mpisize == 1 || forest->ghosts != NULL, - "Ghost structure is needed for t8_forest_leaf_face_neighbors " - "but was not found in forest.\n"); - - if (forest_is_balanced) { - /* In a balanced forest, the leaf neighbor of a leaf is either the neighbor element itself, - * its parent or its children at the face. */ - eclass = t8_forest_get_tree_class (forest, ltreeid); - - if (orientation) { - *orientation = t8_forest_leaf_face_orientation (forest, ltreeid, scheme, leaf, face); - } + T8_ASSERT (pelement_indices != NULL); + T8_ASSERT (dual_faces != NULL); + T8_ASSERT (num_neighbors != NULL); - /* At first we compute these children of the face neighbor elements of leaf. For this, we need the - * neighbor tree's eclass, scheme, and tree id */ - *pneigh_eclass = t8_forest_element_neighbor_eclass (forest, ltreeid, leaf, face); - /* If we are at the maximum refinement level, we compute the neighbor instead */ - at_maxlevel = scheme->element_get_level (eclass, leaf) == t8_forest_get_maxlevel (forest); - if (at_maxlevel) { - num_children_at_face = 1; - neighbor_leaves = *pneighbor_leaves = T8_ALLOC (t8_element_t *, 1); - *dual_faces = T8_ALLOC (int, 1); - scheme->element_new (*pneigh_eclass, num_children_at_face, neighbor_leaves); - /* Compute neighbor element and global treeid of the neighbor */ - gneigh_treeid = t8_forest_element_face_neighbor (forest, ltreeid, leaf, neighbor_leaves[0], *pneigh_eclass, face, - *dual_faces); - } - else { - /* Allocate neighbor element */ - num_children_at_face = scheme->element_get_num_face_children (eclass, leaf, face); - neighbor_leaves = *pneighbor_leaves = T8_ALLOC (t8_element_t *, num_children_at_face); - *dual_faces = T8_ALLOC (int, num_children_at_face); - scheme->element_new (*pneigh_eclass, num_children_at_face, neighbor_leaves); - /* Compute neighbor elements and global treeid of the neighbor */ - gneigh_treeid = t8_forest_element_half_face_neighbors (forest, ltreeid, leaf, neighbor_leaves, *pneigh_eclass, - face, num_children_at_face, *dual_faces); - } - if (gneigh_tree) { - *gneigh_tree = gneigh_treeid; - } - if (gneigh_treeid < 0) { - /* There exists no face neighbor across this face, we return with this info */ - scheme->element_destroy (*pneigh_eclass, num_children_at_face, neighbor_leaves); - T8_FREE (neighbor_leaves); - T8_FREE (*dual_faces); - *dual_faces = NULL; - *num_neighbors = 0; - *pelement_indices = NULL; - *pneighbor_leaves = NULL; - return; - } - T8_ASSERT (gneigh_treeid >= 0 && gneigh_treeid < forest->global_num_trees); - /* We have computed the half face neighbor elements, we now compute their owners, - * if they differ, we know that the half face neighbors are the neighbor leaves. - * If the owners do not differ, we have to check if the neighbor leaf is their - * parent or grandparent. */ - owners = T8_ALLOC (int, num_children_at_face); - different_owners = 0; - have_ghosts = 0; - for (ineigh = 0; ineigh < num_children_at_face; ineigh++) { - /* At first, we check whether the current rank owns the neighbor, since - * this is a constant time check and it is the most common case */ - if (t8_forest_element_check_owner (forest, neighbor_leaves[ineigh], gneigh_treeid, *pneigh_eclass, - forest->mpirank, at_maxlevel)) { - owners[ineigh] = forest->mpirank; - /* The neighbor tree is also a local tree. we store its local treeid */ - lneigh_treeid = t8_forest_get_local_id (forest, gneigh_treeid); - } - else { - owners[ineigh] = t8_forest_element_find_owner (forest, gneigh_treeid, neighbor_leaves[ineigh], *pneigh_eclass); - /* Store that at least one neighbor is a ghost */ - have_ghosts = 1; - } - if (ineigh > 0) { - /* Check if all owners are the same for all neighbors or not */ - different_owners = different_owners || (owners[ineigh] != owners[ineigh - 1]); +#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 (have_ghosts) { - /* At least one neighbor is a ghost, we compute the ghost treeid of the neighbor - * tree. */ - lghost_treeid = t8_forest_ghost_get_ghost_treeid (forest, gneigh_treeid); - T8_ASSERT (lghost_treeid >= 0); + } + + 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); + } + } } - /* TODO: Maybe we do not need to compute the owners. It suffices to know - * whether the neighbor is owned by mpirank or not. */ - - if (!different_owners) { - /* The face neighbors belong to the same process, we thus need to determine - * if they are leaves or their parent or grandparent. */ - neigh_id = scheme->element_get_linear_id (*pneigh_eclass, neighbor_leaves[0], forest->maxlevel); - if (owners[0] != forest->mpirank) { - /* The elements are ghost elements of the same owner */ - const t8_element_array_t *element_array = t8_forest_ghost_get_tree_leaf_elements (forest, lghost_treeid); - /* Find the index in element_array of the leaf ancestor of the first neighbor. - * This is either the neighbor itself or its parent, or its grandparent */ - element_index = t8_forest_bin_search_lower (element_array, neigh_id, forest->maxlevel); - T8_ASSERT (element_index >= 0); - - /* Get the element */ - ancestor = t8_forest_ghost_get_leaf_element (forest, lghost_treeid, element_index); - /* Add the number of ghost elements on previous ghost trees and the number of local elements. */ - element_index += t8_forest_ghost_get_tree_element_offset (forest, lghost_treeid); - element_index += t8_forest_get_local_num_leaf_elements (forest); - T8_ASSERT (forest->local_num_leaf_elements <= element_index - && element_index < forest->local_num_leaf_elements + t8_forest_get_num_ghosts (forest)); + } + + 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 { - /* the elements are local elements */ - const t8_element_array_t *element_array = t8_forest_get_tree_leaf_element_array (forest, lneigh_treeid); - /* Find the index in element_array of the leaf ancestor of the first neighbor. - * This is either the neighbor itself or its parent, or its grandparent */ - element_index = t8_forest_bin_search_lower (element_array, neigh_id, forest->maxlevel); - /* Get the element */ - ancestor = t8_forest_get_tree_leaf_element (t8_forest_get_tree (forest, lneigh_treeid), element_index); - /* Add the element offset of this tree to the index */ - element_index += t8_forest_get_tree_element_offset (forest, lneigh_treeid); - } - if (scheme->element_compare (*pneigh_eclass, ancestor, neighbor_leaves[0]) < 0) { - /* ancestor is a real ancestor, and thus the neighbor is either the parent - * or the grandparent of the half neighbors. We can return it and the indices. */ - /* We need to determine the dual face */ - if (scheme->element_get_level (*pneigh_eclass, ancestor) == scheme->element_get_level (eclass, leaf)) { - /* The ancestor is the same-level neighbor of leaf */ - if (!at_maxlevel) { - /* its dual face is the face of the parent of the first neighbor leaf */ - *dual_faces[0] = scheme->element_face_get_parent_face (*pneigh_eclass, neighbor_leaves[0], *dual_faces[0]); - } + /* 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 { - /* The ancestor is the parent of the parent */ - T8_ASSERT (scheme->element_get_level (*pneigh_eclass, ancestor) - == scheme->element_get_level (eclass, leaf) - 1); - - *dual_faces[0] = scheme->element_face_get_parent_face (*pneigh_eclass, neighbor_leaves[0], *dual_faces[0]); - if (!at_maxlevel) { - /* We need to compute the dual face of the grandparent. */ - /* Construct the parent of the grand child */ - scheme->element_get_parent (*pneigh_eclass, neighbor_leaves[0], neighbor_leaves[0]); - /* Compute the face id of the parent's face */ - *dual_faces[0] = scheme->element_face_get_parent_face (*pneigh_eclass, neighbor_leaves[0], *dual_faces[0]); - } - } - /* free memory */ - scheme->element_destroy (*pneigh_eclass, num_children_at_face - 1, neighbor_leaves + 1); - /* copy the ancestor */ - scheme->element_copy (*pneigh_eclass, ancestor, neighbor_leaves[0]); - /* set return values */ - *num_neighbors = 1; - *pelement_indices = T8_ALLOC (t8_locidx_t, 1); - (*pelement_indices)[0] = element_index; - - T8_FREE (owners); - return; - } - } - /* The leaves are the face neighbors that we are looking for. */ - /* The face neighbors either belong to different processes and thus must be leaves - * in the forest, or the ancestor leaf of the first half neighbor is the half - * neighbor itself and thus all half neighbors must be leaves. - * Since the forest is balanced, we found all neighbor leaves. - * It remains to compute their local ids */ - *num_neighbors = num_children_at_face; - *pelement_indices = T8_ALLOC (t8_locidx_t, num_children_at_face); - element_indices = *pelement_indices; - for (ineigh = 0; ineigh < num_children_at_face; ineigh++) { - /* Compute the linear id at maxlevel of the neighbor leaf */ - neigh_id = scheme->element_get_linear_id (*pneigh_eclass, neighbor_leaves[ineigh], forest->maxlevel); - /* Get a pointer to the element array in which the neighbor lies and search for the element's index in this array. - * This is either the local leaf array of the local tree or the corresponding leaf array in the ghost structure */ - if (owners[ineigh] == forest->mpirank) { - /* The neighbor is a local leaf */ - const t8_element_array_t *element_array = t8_forest_get_tree_leaf_element_array (forest, lneigh_treeid); - /* Find the index of the neighbor in the array */ - element_indices[ineigh] = t8_forest_bin_search_lower (element_array, neigh_id, forest->maxlevel); - T8_ASSERT (element_indices[ineigh] >= 0); - /* We have to add the tree's element offset to the index found to get the actual local element id */ - element_indices[ineigh] += t8_forest_get_tree_element_offset (forest, lneigh_treeid); -#if T8_ENABLE_DEBUG - /* We check whether the element is really the element at this local id */ - { - t8_locidx_t check_ltreeid; - const t8_element_t *check_element - = t8_forest_get_leaf_element (forest, element_indices[ineigh], &check_ltreeid); - T8_ASSERT (check_ltreeid == lneigh_treeid); - T8_ASSERT (scheme->element_is_equal (*pneigh_eclass, check_element, neighbor_leaves[ineigh])); + 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; } -#endif + 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); } - else { - /* The neighbor is a ghost */ - const t8_element_array_t *element_array = t8_forest_ghost_get_tree_leaf_elements (forest, lghost_treeid); - /* Find the index of the neighbor in the array */ - element_indices[ineigh] = t8_forest_bin_search_lower (element_array, neigh_id, forest->maxlevel); - -#if T8_ENABLE_DEBUG - /* We check whether the element is really the element at this local id */ - { - t8_element_t *check_element; - check_element = t8_forest_ghost_get_leaf_element (forest, lghost_treeid, element_indices[ineigh]); - T8_ASSERT (scheme->element_is_equal (*pneigh_eclass, check_element, neighbor_leaves[ineigh])); + // 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]; + } } -#endif - /* Add the element offset of previous ghosts to this index */ - element_indices[ineigh] += t8_forest_ghost_get_tree_element_offset (forest, lghost_treeid); - /* Add the number of all local elements to this index */ - element_indices[ineigh] += t8_forest_get_local_num_leaf_elements (forest); + // 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 for loop over neighbor leaves */ - T8_FREE (owners); + } // End if neighbors exist (first_leaf_index > 0) + // clean up memory allocated with new + delete leaf_array; } - else { - /* TODO: implement unbalanced version */ - SC_ABORT_NOT_REACHED (); + 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 (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *leaf, - t8_element_t **pneighbor_leaves[], int face, int *dual_faces[], int *num_neighbors, - t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass, int forest_is_balanced) +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, forest_is_balanced, NULL, NULL); + 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; - t8_element_t **neighbor_leaves; + 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); @@ -1918,7 +2106,7 @@ t8_forest_print_all_leaf_neighbors (t8_forest_t forest) /* 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, 1); + &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"); @@ -1928,8 +2116,6 @@ t8_forest_print_all_leaf_neighbors (t8_forest_t forest) } t8_debugf ("%s\n", buffer); if (num_neighbors > 0) { - scheme->element_destroy (neigh_eclass, num_neighbors, neighbor_leaves); - T8_FREE (element_indices); T8_FREE (neighbor_leaves); T8_FREE (dual_faces); @@ -2555,7 +2741,12 @@ t8_forest_element_owners_at_neigh_face (t8_forest_t forest, t8_locidx_t ltreeid, int dual_face; t8_gloidx_t neigh_tree; - /* Aallocate memory for the neighbor element */ + if (neigh_class == T8_ECLASS_INVALID) { + /* There is no face neighbor, we indicate this by setting the array to 0 */ + sc_array_resize (owners, 0); + return; + } + /* Allocate memory for the neighbor element */ T8_ASSERT (T8_ECLASS_ZERO <= neigh_class && neigh_class < T8_ECLASS_COUNT); scheme->element_new (neigh_class, 1, &face_neighbor); /* clang-format off */ @@ -2583,6 +2774,12 @@ t8_forest_element_owners_at_neigh_face_bounds (t8_forest_t forest, t8_locidx_t l int dual_face; t8_gloidx_t neigh_tree; + if (neigh_class == T8_ECLASS_INVALID) { + // There is no face neighbor, so there is no owner + *lower = 1; + *upper = 0; + return; + } if (*lower >= *upper) { /* There is no owner or it is unique */ return; diff --git a/src/t8_forest/t8_forest_balance.cxx b/src/t8_forest/t8_forest_balance.cxx index ae67e5dc21..6e8d488a48 100644 --- a/src/t8_forest/t8_forest_balance.cxx +++ b/src/t8_forest/t8_forest_balance.cxx @@ -86,34 +86,39 @@ t8_forest_balance_adapt (t8_forest_t forest, t8_forest_t forest_from, const t8_l for (iface = 0; iface < num_faces; iface++) { /* Get the element class and scheme of the face neighbor */ neigh_class = t8_forest_element_neighbor_eclass (forest_from, ltree_id, element, iface); - /* Allocate memory for the number of half face neighbors */ - num_half_neighbors = scheme->element_get_num_face_children (tree_class, element, iface); - half_neighbors = T8_ALLOC (t8_element_t *, num_half_neighbors); - scheme->element_new (neigh_class, num_half_neighbors, half_neighbors); - /* Compute the half face neighbors of element at this face */ - neighbor_tree = t8_forest_element_half_face_neighbors (forest_from, ltree_id, element, half_neighbors, - neigh_class, iface, num_half_neighbors, NULL); - if (neighbor_tree >= 0) { - /* The face neighbors do exist, check for each one, whether it has + if (neigh_class != T8_ECLASS_INVALID) { + /* Check for each face neighbor, whether it has + * local or ghost leaf descendants in the forest. + * If so, the element will be refined. */ + + /* Allocate memory for the number of half face neighbors */ + num_half_neighbors = scheme->element_get_num_face_children (tree_class, element, iface); + half_neighbors = T8_ALLOC (t8_element_t *, num_half_neighbors); + scheme->element_new (neigh_class, num_half_neighbors, half_neighbors); + /* Compute the half face neighbors of element at this face */ + neighbor_tree = t8_forest_element_half_face_neighbors (forest_from, ltree_id, element, half_neighbors, + neigh_class, iface, num_half_neighbors, NULL); + if (neighbor_tree >= 0) { + /* The face neighbors do exist, check for each one, whether it has * local or ghost leaf descendants in the forest. * If so, the element will be refined. */ - for (ineigh = 0; ineigh < num_half_neighbors; ineigh++) { - if (t8_forest_element_has_leaf_desc (forest_from, neighbor_tree, half_neighbors[ineigh], neigh_class)) { - /* This element should be refined */ - *pdone = 0; - /* clean-up */ - scheme->element_destroy (neigh_class, num_half_neighbors, half_neighbors); - T8_FREE (half_neighbors); - return 1; + for (ineigh = 0; ineigh < num_half_neighbors; ineigh++) { + if (t8_forest_element_has_leaf_desc (forest_from, neighbor_tree, half_neighbors[ineigh], neigh_class)) { + /* This element should be refined */ + *pdone = 0; + /* clean-up */ + scheme->element_destroy (neigh_class, num_half_neighbors, half_neighbors); + T8_FREE (half_neighbors); + return 1; + } } } + /* clean-up */ + scheme->element_destroy (neigh_class, num_half_neighbors, half_neighbors); + T8_FREE (half_neighbors); } - /* clean-up */ - scheme->element_destroy (neigh_class, num_half_neighbors, half_neighbors); - T8_FREE (half_neighbors); } } - return 0; } diff --git a/src/t8_forest/t8_forest_general.h b/src/t8_forest/t8_forest_general.h index 2add53783f..8cde5e01fd 100644 --- a/src/t8_forest/t8_forest_general.h +++ b/src/t8_forest/t8_forest_general.h @@ -3,7 +3,7 @@ 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 + 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 @@ -584,9 +584,9 @@ 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. - * \param [in] forest The forest. Must have a valid ghost layer. - * \param [in] ltreeid A local tree id. +/** 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. @@ -599,17 +599,13 @@ t8_forest_leaf_face_orientation (t8_forest_t forest, const t8_locidx_t ltreeid, * 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 [in] forest_is_balanced True if we know that \a forest is balanced, false - * otherwise. - * \note If there are no face neighbors, then *neighbor_leaves = NULL, num_neighbors = 0, + * \note If there are no face neighbors, then *pneighbor_leaves = NULL, num_neighbors = 0, * and *pelement_indices = NULL on output. - * \note Currently \a forest must be balanced. * \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 along the boundary face. (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) { - * scheme->element_destroy (pneigh_eclass, num_neighbors, pneighbor_leaves); * T8_FREE (pneighbor_leaves); * T8_FREE (pelement_indices); * T8_FREE (dual_faces); @@ -617,14 +613,14 @@ t8_forest_leaf_face_orientation (t8_forest_t forest, const t8_locidx_t ltreeid, * */ void -t8_forest_leaf_face_neighbors (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *leaf, - t8_element_t **pneighbor_leaves[], int face, int *dual_faces[], int *num_neighbors, - t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass, int forest_is_balanced); +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. - * \param [in] leaf A leaf in tree \a ltreeid of \a 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_or_ghost A leaf or ghost leaf element 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 @@ -636,22 +632,18 @@ t8_forest_leaf_face_neighbors (t8_forest_t forest, t8_locidx_t ltreeid, const t8 * 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 [in] forest_is_balanced True if we know that \a forest is balanced, false - * otherwise. * \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 *neighbor_leaves = NULL, num_neighbors = 0, + * \note If there are no face neighbors, then *pneighbor_leaves = NULL, num_neighbors = 0, * and *pelement_indices = NULL on output. - * \note Currently \a forest must be balanced. * \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) { - * scheme->element_destroy (pneigh_eclass, num_neighbors, pneighbor_leaves); * T8_FREE (pneighbor_leaves); * T8_FREE (pelement_indices); * T8_FREE (dual_faces); @@ -659,10 +651,29 @@ t8_forest_leaf_face_neighbors (t8_forest_t forest, t8_locidx_t ltreeid, const t8 * */ void -t8_forest_leaf_face_neighbors_ext (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *leaf, - t8_element_t **pneighbor_leaves[], int face, int *dual_faces[], int *num_neighbors, - t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass, int forest_is_balanced, - t8_gloidx_t *gneigh_tree, int *orientation); +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. @@ -855,17 +866,19 @@ 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 +/** 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 in which the element lies. - * \param [in] elem An element in the tree \a ltreeid. + * \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 local tree id of the tree in which the face - * neighbor of \a elem across \a face lies. + * \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 (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *elem, int face); +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. @@ -882,7 +895,8 @@ t8_forest_element_neighbor_eclass (t8_forest_t forest, t8_locidx_t ltreeid, cons * \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. + * -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, diff --git a/src/t8_forest/t8_forest_ghost.cxx b/src/t8_forest/t8_forest_ghost.cxx index 8acfe43026..15a6d80b89 100644 --- a/src/t8_forest/t8_forest_ghost.cxx +++ b/src/t8_forest/t8_forest_ghost.cxx @@ -864,69 +864,66 @@ t8_forest_ghost_fill_remote (t8_forest_t forest, t8_forest_ghost_t ghost, int gh is_atom = 0; } for (iface = 0; iface < num_faces; iface++) { - /* TODO: Check whether the neighbor element is inside the forest, - * if not then do not compute the half_neighbors. - * This will save computing time. Needs an "element is in forest" function - * Currently we perform this check in the half_neighbors function. */ - /* Get the element class of the neighbor tree */ const t8_eclass_t neigh_class = t8_forest_element_neighbor_eclass (forest, itree, elem, iface); - if (ghost_method == 0) { - /* Use half neighbors */ - /* Get the number of face children of the element at this face */ - num_face_children = scheme->element_get_num_face_children (tree_class, elem, iface); - /* regrow the half_neighbors array if necessary. - * We also need to reallocate it, if the element class of the neighbor - * changes */ - if (max_num_face_children < num_face_children || last_class != neigh_class) { - half_neighbors = T8_ALLOC (t8_element_t *, num_face_children); - /* Allocate memory for the half size face neighbors */ - scheme->element_new (neigh_class, num_face_children, half_neighbors); - max_num_face_children = num_face_children; - last_class = neigh_class; - } - if (!is_atom) { - /* Construct each half size neighbor */ - neighbor_tree = t8_forest_element_half_face_neighbors (forest, itree, elem, half_neighbors, neigh_class, - iface, num_face_children, NULL); - } + if (neigh_class != T8_ECLASS_INVALID) { // Only continue if a face neighbor exists + if (ghost_method == 0) { + /* Use half neighbors */ + /* Get the number of face children of the element at this face */ + num_face_children = scheme->element_get_num_face_children (tree_class, elem, iface); + /* regrow the half_neighbors array if necessary. + * We also need to reallocate it, if the element class of the neighbor + * changes */ + if (max_num_face_children < num_face_children || last_class != neigh_class) { + half_neighbors = T8_ALLOC (t8_element_t *, num_face_children); + /* Allocate memory for the half size face neighbors */ + scheme->element_new (neigh_class, num_face_children, half_neighbors); + max_num_face_children = num_face_children; + last_class = neigh_class; + } + if (!is_atom) { + /* Construct each half size neighbor */ + neighbor_tree = t8_forest_element_half_face_neighbors (forest, itree, elem, half_neighbors, neigh_class, + iface, num_face_children, NULL); + } + else { + int dummy_neigh_face; + /* This element has maximum level, we only construct its neighbor */ + neighbor_tree = t8_forest_element_face_neighbor (forest, itree, elem, half_neighbors[0], neigh_class, + iface, &dummy_neigh_face); + } + if (neighbor_tree >= 0) { + /* If there exist face neighbor elements (we are not at a domain boundary */ + /* Find the owner process of each face_child */ + for (ichild = 0; ichild < num_face_children; ichild++) { + /* find the owner */ + owner = t8_forest_element_find_owner (forest, neighbor_tree, half_neighbors[ichild], neigh_class); + T8_ASSERT (0 <= owner && owner < forest->mpisize); + if (owner != forest->mpirank) { + /* Add the element as a remote element */ + t8_ghost_add_remote (forest, ghost, owner, itree, elem, ielem); + } + } + } + scheme->element_destroy (neigh_class, num_face_children, half_neighbors); + T8_FREE (half_neighbors); + } /* end ghost_method 0 */ else { - int dummy_neigh_face; - /* This element has maximum level, we only construct its neighbor */ - neighbor_tree = t8_forest_element_face_neighbor (forest, itree, elem, half_neighbors[0], neigh_class, iface, - &dummy_neigh_face); - } - if (neighbor_tree >= 0) { - /* If there exist face neighbor elements (we are not at a domain boundary */ - /* Find the owner process of each face_child */ - for (ichild = 0; ichild < num_face_children; ichild++) { - /* find the owner */ - owner = t8_forest_element_find_owner (forest, neighbor_tree, half_neighbors[ichild], neigh_class); + size_t iowner; + /* Construct the owners at the face of the neighbor element */ + t8_forest_element_owners_at_neigh_face (forest, itree, elem, iface, &owners); + /* Iterate over all owners and if any is not the current process, + * add this element as remote */ + for (iowner = 0; iowner < owners.elem_count; iowner++) { + owner = *(int *) sc_array_index (&owners, iowner); T8_ASSERT (0 <= owner && owner < forest->mpisize); if (owner != forest->mpirank) { /* Add the element as a remote element */ t8_ghost_add_remote (forest, ghost, owner, itree, elem, ielem); } } + sc_array_truncate (&owners); } - scheme->element_destroy (neigh_class, num_face_children, half_neighbors); - T8_FREE (half_neighbors); - } /* end ghost_method 0 */ - else { - size_t iowner; - /* Construct the owners at the face of the neighbor element */ - t8_forest_element_owners_at_neigh_face (forest, itree, elem, iface, &owners); - /* Iterate over all owners and if any is not the current process, - * add this element as remote */ - for (iowner = 0; iowner < owners.elem_count; iowner++) { - owner = *(int *) sc_array_index (&owners, iowner); - T8_ASSERT (0 <= owner && owner < forest->mpisize); - if (owner != forest->mpirank) { - /* Add the element as a remote element */ - t8_ghost_add_remote (forest, ghost, owner, itree, elem, ielem); - } - } - sc_array_truncate (&owners); } } /* end face loop */ } /* end element loop */ diff --git a/src/t8_forest/t8_forest_ghost.h b/src/t8_forest/t8_forest_ghost.h index f6936bfbb3..a3897c2060 100644 --- a/src/t8_forest/t8_forest_ghost.h +++ b/src/t8_forest/t8_forest_ghost.h @@ -84,7 +84,7 @@ t8_forest_ghost_tree_num_leaf_elements (t8_forest_t forest, t8_locidx_t lghost_t /** Get a pointer to the ghost leaf element array of a ghost tree. * \param [in] forest The forest. Ghost layer must exist. - * \param [in] lghost_tree The ghost tree id of a ghost tree. + * \param [in] lghost_tree The ghost tree id of a ghost tree. 0 <= \a lghost_tree < num_ghost_trees * \return A pointer to the array of ghost leaf elements of the tree. * \a forest must be committed before calling this function. */ diff --git a/src/t8_forest/t8_forest_iterate.cxx b/src/t8_forest/t8_forest_iterate.cxx index bf90cc9fbf..35a95d51c0 100644 --- a/src/t8_forest/t8_forest_iterate.cxx +++ b/src/t8_forest/t8_forest_iterate.cxx @@ -94,84 +94,103 @@ t8_forest_split_array (const t8_element_t *element, const t8_element_array_t *le void t8_forest_iterate_faces (const t8_forest_t forest, const t8_locidx_t ltreeid, const t8_element_t *element, - const int face, const t8_element_array_t *leaf_elements, void *user_data, - const t8_locidx_t tree_lindex_of_first_leaf, const t8_forest_iterate_face_fn callback) + const int face, const t8_element_array_t *const leaf_elements, + const t8_locidx_t tree_lindex_of_first_leaf, const t8_forest_iterate_face_fn callback, + void *user_data) { - const t8_scheme *scheme = t8_forest_get_scheme (forest); - t8_eclass_t eclass; - t8_element_t **face_children; - int child_face, num_face_children, iface; - int *child_indices; - size_t *split_offsets, indexa, indexb, elem_count; - t8_element_array_t face_child_leaves; - + t8_debugf ("Entering t8_forest_iterate_faces with leaf_index %i and %li total leaves.\n", tree_lindex_of_first_leaf, + t8_element_array_get_count (leaf_elements)); T8_ASSERT (t8_forest_is_committed (forest)); - T8_ASSERT (0 <= ltreeid && ltreeid < t8_forest_get_num_local_trees (forest)); + const t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (forest); +#if T8_ENABLE_DEBUG + const t8_locidx_t num_ghost_trees = t8_forest_get_num_ghost_trees (forest); + const t8_locidx_t num_local_and_ghost_trees = num_local_trees + num_ghost_trees; + T8_ASSERT (0 <= ltreeid && ltreeid < num_local_and_ghost_trees); +#endif + + // Check whether we are in a local tree or ghost tree + const bool tree_is_local = t8_forest_tree_is_local (forest, ltreeid); + const bool tree_is_ghost = !tree_is_local; + const t8_locidx_t local_or_ghost_tree_id = tree_is_local ? ltreeid : ltreeid - num_local_trees; - elem_count = t8_element_array_get_count (leaf_elements); + const size_t elem_count = t8_element_array_get_count (leaf_elements); if (elem_count == 0) { /* There are no leaves left, so we have nothing to do */ return; } - eclass = t8_forest_get_tree_class (forest, ltreeid); + const t8_eclass_t eclass = t8_forest_get_tree_class (forest, ltreeid); + const t8_scheme *scheme = t8_forest_get_scheme (forest); - if (elem_count == 1) { - /* There is only one leaf left, we check whether it is the same as element - * and if so call the callback function */ + // TODO: This does a costly binary search. Is there a better criterion to check + // whether element is a leaf or not? + // We tried (elem_count == 1) but this does not work since we could + // start the search with a full family and element being one sibling. + // In that case, element is a leaf but elem_count is not 1. + bool is_leaf = t8_forest_element_is_leaf_or_ghost (forest, element, local_or_ghost_tree_id, tree_is_ghost); + +#if T8_ENABLE_DEBUG + if (!is_leaf) { + /* Check whether element has smaller level than the first leaf */ const t8_element_t *leaf = t8_element_array_index_locidx (leaf_elements, 0); - T8_ASSERT (t8_forest_element_is_leaf (forest, leaf, ltreeid)); - if (scheme->element_is_equal (eclass, element, leaf)) { - /* The element is the leaf, we are at the last stage of the recursion - * and can call the callback. */ - (void) callback (forest, ltreeid, leaf, face, user_data, tree_lindex_of_first_leaf); - return; + T8_ASSERT (t8_forest_element_is_leaf_or_ghost (forest, leaf, local_or_ghost_tree_id, tree_is_ghost)); + T8_ASSERT (scheme->element_get_level (eclass, element) < scheme->element_get_level (eclass, leaf)); + + // Verify that all leaves in leaf_elements are descendants of element + for (t8_locidx_t ileaf_index = 0; (size_t) ileaf_index < elem_count; ++ileaf_index) { + const t8_element_t *ileaf = t8_element_array_index_locidx (leaf_elements, ileaf_index); + T8_ASSERT (scheme->element_is_ancestor (eclass, element, ileaf)); } } -#if T8_ENABLE_DEBUG - /* Check whether element has greater level than the first leaf */ - const t8_element_t *leaf = t8_element_array_index_locidx (leaf_elements, 0); - T8_ASSERT (t8_forest_element_is_leaf (forest, leaf, ltreeid)); - T8_ASSERT (scheme->element_get_level (eclass, element) < scheme->element_get_level (eclass, leaf)); #endif - /* Call the callback function element, we pass -index - 1 as index to indicate - * element is not a leaf, if it returns true, we continue with the top-down recursion */ - if (callback (forest, ltreeid, element, face, user_data, -tree_lindex_of_first_leaf - 1)) { - /* Enter the recursion */ - /* We compute all face children of E, compute their leaf arrays and call iterate_faces */ - /* allocate the memory to store the face children */ - num_face_children = scheme->element_get_num_face_children (eclass, element, face); - face_children = T8_ALLOC (t8_element_t *, num_face_children); - scheme->element_new (eclass, num_face_children, face_children); - /* Memory for the child indices of the face children */ - child_indices = T8_ALLOC (int, num_face_children); - /* Memory for the indices that split the leaf_elements array */ - split_offsets = T8_ALLOC (size_t, scheme->element_get_num_children (eclass, element) + 1); - /* Compute the face children */ - scheme->element_get_children_at_face (eclass, element, face, face_children, num_face_children, child_indices); - /* Split the leaves array in portions belonging to the children of element */ - t8_forest_split_array (element, leaf_elements, split_offsets); - for (iface = 0; iface < num_face_children; iface++) { - /* Check if there are any leaf elements for this face child */ - indexa = split_offsets[child_indices[iface]]; /* first leaf of this face child */ - indexb = split_offsets[child_indices[iface] + 1]; /* first leaf of next child */ - if (indexa < indexb) { - /* There exist leaves of this face child in leaf_elements, - * we construct an array of these leaves */ - t8_element_array_init_view (&face_child_leaves, leaf_elements, indexa, indexb - indexa); - /* Compute the corresponding face number of this face child */ - child_face = scheme->element_face_get_child_face (eclass, element, face, iface); - /* Enter the recursion */ - t8_forest_iterate_faces (forest, ltreeid, face_children[iface], child_face, &face_child_leaves, user_data, - indexa + tree_lindex_of_first_leaf, callback); - } + /* Call the callback function element, if it returns true, we continue with the top-down recursion */ + const int ret + = callback (forest, ltreeid, element, face, is_leaf, leaf_elements, tree_lindex_of_first_leaf, user_data); + if (!ret || is_leaf) { + // The callback returned false or the element is a leaf. + // We abort the recursion. + return; + } + + /* Enter the recursion */ + /* We compute all face children of E, compute their leaf arrays and call iterate_faces */ + /* allocate the memory to store the face children */ + const int num_face_children = scheme->element_get_num_face_children (eclass, element, face); + t8_element_t **face_children = T8_ALLOC (t8_element_t *, num_face_children); + scheme->element_new (eclass, num_face_children, face_children); + /* Memory for the child indices of the face children */ + int *child_indices = T8_ALLOC (int, num_face_children); + /* Memory for the indices that split the leaf_elements array */ + const int num_children = scheme->element_get_num_children (eclass, element); + size_t *split_offsets = T8_ALLOC (size_t, num_children + 1); + /* Compute the face children */ + scheme->element_get_children_at_face (eclass, element, face, face_children, num_face_children, child_indices); + /* Split the leaves array in portions belonging to the children of element */ + t8_forest_split_array (element, leaf_elements, split_offsets); + for (int iface = 0; iface < num_face_children; iface++) { + /* Check if there are any leaf elements for this face child */ + T8_ASSERT (0 <= child_indices[iface]); + T8_ASSERT (child_indices[iface] < num_children + 1); + const size_t indexa = split_offsets[child_indices[iface]]; /* first leaf of this face child */ + const size_t indexb = split_offsets[child_indices[iface] + 1]; /* first leaf of next child */ + t8_debugf ("Computed indices for face child %i: %zd %zd\n", iface, indexa, indexb); + if (indexa < indexb) { + /* There exist leaves of this face child in leaf_elements, + * we construct an array of these leaves */ + t8_element_array_t face_child_leaves; + t8_element_array_init_view (&face_child_leaves, leaf_elements, indexa, indexb - indexa); + /* Compute the corresponding face number of this face child */ + const int child_face = scheme->element_face_get_child_face (eclass, element, face, iface); + /* Enter the recursion */ + t8_forest_iterate_faces (forest, ltreeid, face_children[iface], child_face, &face_child_leaves, + indexa + tree_lindex_of_first_leaf, callback, user_data); } - /* clean-up */ - scheme->element_destroy (eclass, num_face_children, face_children); - T8_FREE (face_children); - T8_FREE (child_indices); - T8_FREE (split_offsets); } + /* clean-up */ + scheme->element_destroy (eclass, num_face_children, face_children); + T8_FREE (face_children); + T8_FREE (child_indices); + T8_FREE (split_offsets); } /** diff --git a/src/t8_forest/t8_forest_iterate.h b/src/t8_forest/t8_forest_iterate.h index 548fb4cf99..06a92e6a08 100644 --- a/src/t8_forest/t8_forest_iterate.h +++ b/src/t8_forest/t8_forest_iterate.h @@ -38,14 +38,17 @@ * \param[in] ltreeid Local index of the tree containing the \a element. * \param[in] element The considered element. * \param[in] face The integer index of the considered face of \a element. + * \param[in] is_leaf True if and only if the currently considered element is a leaf element. + * \param[in] leaf_elements The array of leaf elements that are descendants of \a element. Sorted by linear index. + * \param[in] tree_leaf_index Tree-local index of the first leaf. * \param[in] user_data Some user-defined data, as void pointer. - * \param[in] tree_leaf_index Tree-local index the first leaf. * * \return Nonzero if the element may touch the face and the top-down search shall be continued, zero otherwise. */ typedef int (*t8_forest_iterate_face_fn) (const t8_forest_t forest, const t8_locidx_t ltreeid, - const t8_element_t *element, const int face, void *user_data, - const t8_locidx_t tree_leaf_index); + const t8_element_t *element, const int face, const int is_leaf, + const t8_element_array_t *leaf_elements, const t8_locidx_t tree_leaf_index, + void *user_data); /** * A call-back function used by \ref t8_forest_search describing a search-criterion. Is called on an element and the @@ -150,6 +153,10 @@ T8_EXTERN_C_BEGIN (); void t8_forest_split_array (const t8_element_t *element, const t8_element_array_t *leaf_elements, size_t *offsets); +/* TODO: comment */ +// TODO: Test this function. Uniform mesh test. Refine always same corner, know that neighbors follow geometric series. +// TODO: user data should be a template parameter in the long run +// TODO: adapt to search interface /** * Iterate over all leaves of an element that touch a given face of the element. * Callback is called in each recursive step with element as input. @@ -164,14 +171,17 @@ t8_forest_split_array (const t8_element_t *element, const t8_element_array_t *le * \param[in] element The considered element. * \param[in] face The integer index of the considered face of \a element. * \param[in] leaf_elements The array of leaf elements that are descendants of \a element. Sorted by linear index. - * \param[in] user_data The user data passed to the \a callback function. - * \param[in] tree_lindex_of_first_leaf Tree-local index of the first leaf. + * \param[in] tree_lindex_of_first_leaf Index of the first leaf of \a element in the tree's leaves. + * The corresponding leaf does not necessarily lie on the face of \a element. + * \note \a tree_lindex_of_first_leaf is not an index in \a leaf_elements. \a leaf_elements may only be a part of the tree's leaves. * \param[in] callback The callback function. + * \param[in] user_data The user data passed to the \a callback function. */ void t8_forest_iterate_faces (const t8_forest_t forest, const t8_locidx_t ltreeid, const t8_element_t *element, - const int face, const t8_element_array_t *leaf_elements, void *user_data, - const t8_locidx_t tree_lindex_of_first_leaf, const t8_forest_iterate_face_fn callback); + const int face, const t8_element_array_t *const leaf_elements, + const t8_locidx_t tree_lindex_of_first_leaf, const t8_forest_iterate_face_fn callback, + void *user_data); /** * Perform a top-down search of the forest, executing a callback on each diff --git a/src/t8_forest/t8_forest_private.h b/src/t8_forest/t8_forest_private.h index 0414a343e8..46dbedb218 100644 --- a/src/t8_forest/t8_forest_private.h +++ b/src/t8_forest/t8_forest_private.h @@ -362,7 +362,7 @@ t8_forest_element_owners_bounds (t8_forest_t forest, t8_gloidx_t gtreeid, const * on output a (better) bound. * * \note If on input \a lower >= \a upper, then the bounds are not changed by this - * algorithm. We interpret \a lower = \a such that the owner is unique and equals \a lower. + * algorithm. We interpret \a lower = \a upper such that the owner is unique and equals \a lower. * \note \a forest must be committed before calling this function. */ void diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 05301f0dcd..aeb71b770e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -156,6 +156,7 @@ add_t8_cpp_test( NAME t8_gtest_ghost_and_owner_parallel SOURCES t8_for add_t8_cpp_test( NAME t8_gtest_balance_parallel SOURCES t8_forest/t8_gtest_balance.cxx ) add_t8_cpp_test( NAME t8_gtest_forest_commit_parallel SOURCES t8_forest/t8_gtest_forest_commit.cxx ) add_t8_cpp_test( NAME t8_gtest_forest_face_normal_serial SOURCES t8_forest/t8_gtest_forest_face_normal.cxx ) +add_t8_cpp_test( NAME t8_gtest_forest_leaf_face_neighbor_parallel SOURCES t8_forest/t8_gtest_leaf_face_neighbors.cxx t8_gtest_adapt_callbacks.cxx ) add_t8_cpp_test( NAME t8_gtest_partition_data_parallel SOURCES t8_forest/t8_gtest_partition_data.cxx ) add_t8_cpp_test( NAME t8_gtest_set_partition_offset_parallel SOURCES t8_forest/t8_gtest_set_partition_offset.cxx ) add_t8_cpp_test( NAME t8_gtest_partition_for_coarsening_parallel SOURCES t8_forest/t8_gtest_partition_for_coarsening.cxx ) diff --git a/test/mesh_handle/t8_gtest_ghost.cxx b/test/mesh_handle/t8_gtest_ghost.cxx index ef028f8fed..4b602b77a3 100644 --- a/test/mesh_handle/t8_gtest_ghost.cxx +++ b/test/mesh_handle/t8_gtest_ghost.cxx @@ -132,14 +132,13 @@ TEST_P (t8_mesh_ghost_test, compare_neighbors_to_forest) EXPECT_EQ (mesh_iterator->get_num_faces (), num_faces); for (int iface = 0; iface < num_faces; iface++) { // --- Get neighbors from forest. --- - t8_element_t** neighbors; + const t8_element_t** neighbors; int num_neighbors; - const int forest_is_balanced = t8_forest_is_balanced (forest); t8_eclass_t neigh_eclass; int* dual_faces; t8_locidx_t* neigh_ids; t8_forest_leaf_face_neighbors (forest, itree, elem, &neighbors, iface, &dual_faces, &num_neighbors, &neigh_ids, - &neigh_eclass, forest_is_balanced); + &neigh_eclass); // --- Get neighbors from mesh element. --- std::vector dual_faces_handle; auto neighbors_handle = mesh_iterator->get_face_neighbors (iface, dual_faces_handle); @@ -155,7 +154,6 @@ TEST_P (t8_mesh_ghost_test, compare_neighbors_to_forest) } // Free memory. if (num_neighbors > 0) { - scheme->element_destroy (neigh_eclass, num_neighbors, neighbors); T8_FREE (neigh_ids); T8_FREE (neighbors); T8_FREE (dual_faces); diff --git a/test/t8_forest/t8_gtest_element_is_leaf.cxx b/test/t8_forest/t8_gtest_element_is_leaf.cxx index d162707a42..c297f2f856 100644 --- a/test/t8_forest/t8_gtest_element_is_leaf.cxx +++ b/test/t8_forest/t8_gtest_element_is_leaf.cxx @@ -30,6 +30,7 @@ #include "test/t8_cmesh_generator/t8_cmesh_example_sets.hxx" #include #include +#include /* In this test we check the t8_forest_element_is_leaf function. * Iterating over all cmesh test cases, we creat a uniform and an adaptive forest. @@ -45,7 +46,8 @@ #define T8_IS_LEAF_MAX_LVL 3 #endif -class element_is_leaf_or_ghost: public testing::TestWithParam> { +struct element_is_leaf_or_ghost: public testing::TestWithParam> +{ protected: void SetUp () override @@ -64,7 +66,6 @@ class element_is_leaf_or_ghost: public testing::TestWithParam 0) { - scheme->element_destroy (neigh_eclass, num_neighbors, neighbors); T8_FREE (neigh_ids); T8_FREE (neighbors); T8_FREE (dual_faces); diff --git a/test/t8_forest/t8_gtest_half_neighbors.cxx b/test/t8_forest/t8_gtest_half_neighbors.cxx index b40fe34212..553a4c90cf 100644 --- a/test/t8_forest/t8_gtest_half_neighbors.cxx +++ b/test/t8_forest/t8_gtest_half_neighbors.cxx @@ -74,19 +74,18 @@ TEST_P (forest_half_neighbors, test_half_neighbors) for (int face = 0; face < scheme->element_get_num_faces (eclass, element); face++) { /* Get the eclass of the face neighbor and get the scheme */ const t8_eclass_t neigh_class = t8_forest_element_neighbor_eclass (forest, itree, element, face); - /* allocate memory for element's neighbor and construct it */ - scheme->element_new (neigh_class, 1, &neighbor); - const t8_locidx_t neigh_tree - = t8_forest_element_face_neighbor (forest, itree, element, neighbor, neigh_class, face, &dual_face); - if (neigh_tree >= 0) { + if (neigh_class != T8_ECLASS_INVALID) { + /* allocate memory for element's neighbor and construct it */ + scheme->element_new (neigh_class, 1, &neighbor); + (void) t8_forest_element_face_neighbor (forest, itree, element, neighbor, neigh_class, face, &dual_face); const int num_face_neighs = scheme->element_get_num_face_children (neigh_class, neighbor, dual_face); t8_element_t **half_neighbors = T8_TESTSUITE_ALLOC (t8_element_t *, num_face_neighs); scheme->element_new (neigh_class, num_face_neighs, half_neighbors); t8_forest_element_half_face_neighbors (forest, itree, element, half_neighbors, neigh_class, face, num_face_neighs, NULL); - /* We now check whether the face children of neighbor are the half neighbors. */ T8_ASSERT (num_face_neighs == scheme->element_get_num_face_children (neigh_class, neighbor, dual_face)); + EXPECT_NE (neigh_class, T8_ECLASS_INVALID); t8_element_t **neighbor_face_children = T8_TESTSUITE_ALLOC (t8_element_t *, num_face_neighs); scheme->element_new (neigh_class, num_face_neighs, neighbor_face_children); int *child_ids = T8_TESTSUITE_ALLOC (int, num_face_neighs); @@ -98,12 +97,12 @@ TEST_P (forest_half_neighbors, test_half_neighbors) << "ineigh = " << ineigh << " face = " << face; } scheme->element_destroy (neigh_class, num_face_neighs, neighbor_face_children); - scheme->element_destroy (neigh_class, num_face_neighs, half_neighbors); - T8_TESTSUITE_FREE (child_ids); T8_TESTSUITE_FREE (neighbor_face_children); + T8_TESTSUITE_FREE (child_ids); + scheme->element_destroy (neigh_class, 1, &neighbor); + scheme->element_destroy (neigh_class, num_face_neighs, half_neighbors); T8_TESTSUITE_FREE (half_neighbors); } - scheme->element_destroy (neigh_class, 1, &neighbor); } } } diff --git a/test/t8_forest/t8_gtest_leaf_face_neighbors.cxx b/test/t8_forest/t8_gtest_leaf_face_neighbors.cxx new file mode 100644 index 0000000000..d1c8e43b95 --- /dev/null +++ b/test/t8_forest/t8_gtest_leaf_face_neighbors.cxx @@ -0,0 +1,574 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2026 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "test/t8_cmesh_generator/t8_cmesh_example_sets.hxx" +#include "t8_test_data_dir.h" +#include + +bool +test_face_neighbors_skip_cmesh (const t8_cmesh_t cmesh) +{ + // We skip empty cmeshes. + return t8_cmesh_is_empty (cmesh); +} + +class forest_face_neighbors: public testing::TestWithParam > { + protected: + void + SetUp () override + { + const int scheme_id = std::get<0> (GetParam ()); + t8_cmesh_t cmesh = std::get<1> (GetParam ())->cmesh_create (); + if (test_face_neighbors_skip_cmesh (cmesh)) { + /* we skip empty cmeshes case */ + t8_cmesh_unref (&cmesh); + GTEST_SKIP (); + } + const t8_scheme *scheme = create_from_scheme_id (scheme_id); + const int uniform_level = 0; + const int adapt_levels = 2; + const int max_adapt_level = uniform_level + adapt_levels; + const bool do_ghost = true; + const bool do_recursive_adapt = true; + forests[0] = t8_forest_new_uniform (cmesh, scheme, uniform_level, do_ghost, sc_MPI_COMM_WORLD); + cmesh = t8_forest_get_cmesh (forests[0]); + t8_forest_ref (forests[0]); + forests[1] = t8_forest_new_adapt (forests[0], t8_test_adapt_first_child, do_recursive_adapt, do_ghost, + (void *) &max_adapt_level); + t8_forest_ref (forests[0]); + // Add another adapted forest that does not create a ghost layer to test + // the face neighbor algorithm on forests without ghost. + const bool dont_do_ghost = false; + forests[2] = t8_forest_new_adapt (forests[0], t8_test_adapt_first_child, do_recursive_adapt, dont_do_ghost, + (void *) &max_adapt_level); + } + + void + TearDown () override + { + for (auto &forest : forests) { + if (forest != nullptr) { + t8_forest_unref (&forest); + } + } + } + + t8_forest_t forests[3] { nullptr, nullptr, nullptr }; +}; + +/* Check that a leaf/ghost element index matches a given leaf/ghost element. + * Returns true if the element with index "element_index" (0<= element_index < num_leaves + num_ghosts) + * lies in the tree "gtreeid" and matches the element "element". */ +static void +verify_leaf_element_index (const t8_forest_t forest, const t8_gloidx_t gtreeid, const t8_locidx_t element_index, + const t8_element_t *element) +{ + const t8_locidx_t num_local_elements = t8_forest_get_local_num_leaf_elements (forest); + + if (element_index < num_local_elements) { + // The element index belongs to a local leaf element. + // Check that the tree is a local tree. + const t8_locidx_t ltreeid = t8_forest_get_local_id (forest, gtreeid); + EXPECT_GE (ltreeid, 0) << "Tree of local element is not a local tree."; + t8_locidx_t check_local_tree_id; + const t8_element_t *element_from_index = t8_forest_get_leaf_element (forest, element_index, &check_local_tree_id); + EXPECT_EQ (check_local_tree_id, ltreeid) << "Element index does not match local tree index."; + EXPECT_EQ (element_from_index, element) + << "Element at index " << element_index << " does not match the given element."; + } + else { + // The element index belongs to a ghost leaf element. + // Check that the neighbor tree is a ghost tree. + const t8_locidx_t ghost_tree_id = t8_forest_ghost_get_ghost_treeid (forest, gtreeid); + EXPECT_GE (ghost_tree_id, 0) << "Tree of ghost element is not a ghost tree."; + // There is no get_ghost_leaf_element function that takes only the element index. + // So we have to convert the element index into a tree local index first. + // Subtract the number of local elements and the number of ghosts in previous trees. + const t8_locidx_t ghost_in_tree_index + = element_index - num_local_elements - t8_forest_ghost_get_tree_element_offset (forest, ghost_tree_id); + const t8_element_t *ghost_element_from_index + = t8_forest_ghost_get_leaf_element (forest, ghost_tree_id, ghost_in_tree_index); + EXPECT_EQ (ghost_element_from_index, element) + << "Neighbor neighbor ghost element at index " << element_index << " is not original element."; + } +} + +TEST_P (forest_face_neighbors, test_face_neighbors) +{ + /* iterate over all elements */ + bool forest_is_uniform = true; // The first forest is uniform. We set this to false at the end of the for loop. + for (auto &forest : forests) { + const t8_cmesh_t cmesh = t8_forest_get_cmesh (forest); + const bool has_ghost = forest->ghosts != NULL; +#if T8_ENABLE_DEBUG + if (t8_cmesh_get_tree_geometry (cmesh, 0) != NULL) { + // Debug vtk output, only if cmesh has a registered geometry + t8_forest_write_vtk (forest, "debug_face_neigh"); + t8_debugf ("writing forest to \'debug_face_neigh\'"); + } +#endif + const t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (forest); + const t8_locidx_t num_ghost_trees = t8_forest_get_num_ghost_trees (forest); + const t8_locidx_t num_local_elements = t8_forest_get_local_num_leaf_elements (forest); + t8_locidx_t ielement_index = 0; + for (t8_locidx_t itree = 0; itree < num_local_trees + num_ghost_trees; itree++) { + const t8_gloidx_t gtree_id = t8_forest_global_tree_id (forest, itree); + const bool is_ghost = itree >= num_local_trees; + const t8_locidx_t ghost_tree_id = itree - num_local_trees; + /* Get the leaf element array */ + const t8_element_array_t *leaf_elements = is_ghost + ? t8_forest_ghost_get_tree_leaf_elements (forest, ghost_tree_id) + : t8_forest_get_tree_leaf_element_array (forest, itree); + const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, itree); + const t8_scheme *scheme = t8_forest_get_scheme (forest); + const t8_locidx_t num_leaves = t8_element_array_get_count (leaf_elements); + const t8_locidx_t cmesh_tree = t8_forest_ltreeid_to_cmesh_ltreeid (forest, itree); + for (t8_locidx_t ileaf = 0; ileaf < num_leaves; ++ileaf, ++ielement_index) { + // Iterate over each leaf element + const t8_element_t *element = t8_element_array_index_locidx (leaf_elements, ileaf); + const int num_faces = scheme->element_get_num_faces (tree_class, element); + for (int iface = 0; iface < num_faces; ++iface) { + // Iterate over all faces and compute the face neighbors + + // preparation + const t8_element_t **neighbor_leaves; + int *dual_faces; + int num_neighbors = 0; + t8_locidx_t *element_indices; + t8_eclass_t neigh_class; + t8_gloidx_t gneigh_tree; + int orientation; + + t8_debugf ("Compute face neighbor for tree %i (%s) element %i (index %i), at face %i.\n", itree, + is_ghost ? "ghost" : "local", ileaf, ielement_index, iface); + + // Actual computation of the face neighbors + t8_forest_leaf_face_neighbors_ext (forest, itree, element, &neighbor_leaves, iface, &dual_faces, + &num_neighbors, &element_indices, &neigh_class, &gneigh_tree, + &orientation); + + t8_debugf ("Tree %i element %i at face %i has %i face neighbors.\n", itree, ileaf, iface, num_neighbors); + + if (gneigh_tree < 0) { + // If there is no neighbor tree then there cannot be any face neighbors. + // Note that there can also be no face neighbors computed if a neighbor tree exists, but + // the element is a ghost and the neighbor element is neither a local element nor ghost. + ASSERT_EQ (num_neighbors, 0); + } + if (num_neighbors == 0) { + // No neighbors are found, check for correctly set return values + ASSERT_TRUE (element_indices == NULL); + ASSERT_TRUE (neighbor_leaves == NULL); + ASSERT_TRUE (dual_faces == NULL); + } + else { + ASSERT_GE (num_neighbors, 0); + ASSERT_TRUE (neighbor_leaves != NULL); + ASSERT_TRUE (element_indices != NULL); + ASSERT_TRUE (dual_faces != NULL); + } + + // Checking for: + // uniform and adapted forest: + // - inner local element has >= 1 face neighbors (= 1 for uniform) + // - inner ghost element has 0 or 1 face neighbors (= 1 for uniform) + // - boundary element has 0 face neighbors + // - If E face f has neighbor E' face f', then + // E' face f' must have neighbor E face f. + + // Now checking for inner and boundary elements. + + // Compute whether this element is a boundary element or not. + // An element is a boundary element if it lies on the tree boundary + // and if the corresponding tree face is at the domain boundary. + // TODO: Use t8_forest_leaf_is_boundary when https://github.com/DLR-AMR/t8code/pull/1081 is integrated + const bool is_root_boundary = scheme->element_is_root_boundary (tree_class, element, iface); + const int tree_face = scheme->element_get_tree_face (tree_class, element, iface); + const bool is_boundary_element + = is_root_boundary && t8_cmesh_tree_face_is_boundary (cmesh, cmesh_tree, tree_face); + + if (!is_boundary_element) { + if (!is_ghost) { // Local element + if (forest_is_uniform) { + // In a uniform forest we must have exactly 1 neighbor. + EXPECT_EQ (num_neighbors, 1) + << "Inner local element should have exactly 1 neighbor, has " << num_neighbors << "."; + } + else if (has_ghost) { + // In an adaptive forest we have 1 or more neighbors. + EXPECT_GE (num_neighbors, 1) + << "Inner local element should have at least 1 neighbor, has " << num_neighbors << "."; + } + } + else { // Ghost element + if (forest_is_uniform) { + // In a uniform forest a ghost element has none or one neighbor. + EXPECT_TRUE (num_neighbors == 0 || num_neighbors == 1) + << "Inner ghost element should have exactly 1 or 0 neighbors, has " << num_neighbors << "."; + } + else if (has_ghost) { + // In an adaptive forest a ghost element has 0 or more neighbors. + EXPECT_GE (num_neighbors, 0) + << "Inner ghost element should have 0 or more neighbors, has " << num_neighbors << "."; + } + } + } + else { + EXPECT_EQ (num_neighbors, 0) << "Boundary element should have exactly 0 neighbors, has " << num_neighbors + << "."; + } + + if (forest_is_uniform) { + ASSERT_TRUE (num_neighbors == 0 || num_neighbors == 1); + // Check the index computation function and that it computes the correct neighbor index. + int check_dual_face; + const t8_locidx_t check_same_level_index = t8_forest_same_level_leaf_face_neighbor_index ( + forest, ielement_index, iface, gtree_id, &check_dual_face); + + if (check_dual_face < 0) { + EXPECT_EQ (num_neighbors, 0); + } + if (check_dual_face >= 0) { + EXPECT_EQ (dual_faces[0], check_dual_face); + EXPECT_EQ (element_indices[0], check_same_level_index); + } + } + + // Check that the neighbor of the neighbor is the original element. + for (int ineigh = 0; ineigh < num_neighbors; ++ineigh) { + const t8_element_t *neighbor = neighbor_leaves[ineigh]; + const int dual_face = dual_faces[ineigh]; + const t8_locidx_t neigh_index = element_indices[ineigh]; + + t8_debugf ("Checking neighbor element %p in (global) tree %li.\n", (void *) neighbor, gneigh_tree); + t8_debugf ("dual face is %i, index is %i\n", dual_face, neigh_index); + +#if T8_ENABLE_DEBUG + scheme->element_debug_print (neigh_class, neighbor); + ASSERT_TRUE (scheme->element_is_valid (neigh_class, neighbor)) + << "Neighbor element " << ineigh << " is not valid"; +#endif + // Check that neighbor index correctly yields neighbor element. + verify_leaf_element_index (forest, gneigh_tree, neigh_index, neighbor); + + // Compute the local treeid of the neighbor tree. + const t8_locidx_t neigh_ltreeid + = neigh_index < num_local_elements + ? gneigh_tree - t8_forest_get_first_local_tree_id (forest) + : t8_forest_ghost_get_ghost_treeid (forest, gneigh_tree) + num_local_trees; + + // preparation + const t8_element_t **neigh_neighbor_leaves; + int *neigh_dual_faces; + int neigh_num_neighbors = 0; + t8_locidx_t *neigh_element_indices; + t8_eclass_t neigh_neigh_class; + t8_gloidx_t neigh_gneigh_tree; + int neigh_orientation; + // Actual computation of the neighbor's face neighbors + t8_forest_leaf_face_neighbors_ext (forest, neigh_ltreeid, neighbor, &neigh_neighbor_leaves, dual_face, + &neigh_dual_faces, &neigh_num_neighbors, &neigh_element_indices, + &neigh_neigh_class, &neigh_gneigh_tree, &neigh_orientation); + +#if T8_ENABLE_DEBUG + t8_debugf ("original element\n"); + scheme->element_debug_print (tree_class, element); + t8_debugf ("neighbor element\n"); + scheme->element_debug_print (neigh_class, neighbor); +#endif + // We must have found at least one face neighbor, namely the original element. + EXPECT_GE (neigh_num_neighbors, 1); + // The neighbor's neighbor tree must be the current tree + EXPECT_EQ (gtree_id, neigh_gneigh_tree); + // The neighbor's scheme must be the current scheme + EXPECT_EQ (tree_class, neigh_neigh_class); + // The neighbor's orientation must be the orientation + EXPECT_EQ (orientation, neigh_orientation); + + // We now (try to) find the original element among the neighbors. + // If it does not exist there was an error. + // If it exists we check that dual face and index were computed correctly. + + int position_of_original_element = -1; + bool found_original = false; + t8_debugf ("Checking all %i neighbors of neighbor for original element:\n", neigh_num_neighbors); + for (int ineighneigh = 0; ineighneigh < neigh_num_neighbors && !found_original; ++ineighneigh) { + // Check that the neighbor of the neighbor element is the original element + const t8_element_t *neigh_of_neigh = neigh_neighbor_leaves[ineighneigh]; + t8_debugf ("Checking neighbor of neighbor #%i:\n", ineighneigh); +#if T8_ENABLE_DEBUG + scheme->element_debug_print (tree_class, neigh_of_neigh); +#endif + if (scheme->element_is_equal (tree_class, element, neigh_of_neigh)) { + position_of_original_element = ineighneigh; + found_original = true; // Stop the for loop + } + } + // We must have found the original element among the neighbors. + ASSERT_TRUE (found_original) << "The original element was not a neighbor of its neighbor."; + + // Check that the dual face of the dual face is the original face + const int neigh_dual_face = neigh_dual_faces[position_of_original_element]; + EXPECT_EQ (neigh_dual_face, iface); + + // Check that the index is correct, i.e. when getting the neighbor neighbor element from the index + // we retrieve the original element. + const t8_locidx_t element_index = neigh_element_indices[position_of_original_element]; + EXPECT_GE (element_index, 0); + + verify_leaf_element_index (forest, neigh_gneigh_tree, element_index, element); + + // Devnote: We are not using T8_TESTSUITE_FREE here, since the memory is allocated inside of t8_forest_leaf_face_neighbors_ext + // by t8code itself using T8_ALLOC. Thus, we need to call T8_FREE to free it. + // clean-up neighbor's neighbors + if (neigh_num_neighbors > 0) { + T8_FREE (neigh_neighbor_leaves); + T8_FREE (neigh_element_indices); + T8_FREE (neigh_dual_faces); + } + } + + // clean-up original element neighbors + if (num_neighbors > 0) { + T8_FREE (neighbor_leaves); + T8_FREE (element_indices); + T8_FREE (dual_faces); + } + } + } + } + forest_is_uniform = false; + } +} + +INSTANTIATE_TEST_SUITE_P (t8_gtest_face_neighbors, forest_face_neighbors, + testing::Combine (AllSchemeCollections, AllCmeshsParam), pretty_print_base_example_scheme); + +/** Adapt callback that refines all elements in the tree with global id 0. + * This callback belong to the test case forest_face_neighbors_two_quad_mesh. + */ +int +t8_test_adapt_first_tree (t8_forest_t forest, [[maybe_unused]] t8_forest_t forest_from, + [[maybe_unused]] t8_locidx_t which_tree, const t8_eclass_t eclass, + [[maybe_unused]] t8_locidx_t lelement_id, const t8_scheme *scheme, + [[maybe_unused]] const int is_family, [[maybe_unused]] const int num_elements, + t8_element_t *elements[]) +{ + T8_ASSERT (!is_family || (is_family && num_elements == scheme->element_get_num_children (eclass, elements[0]))); + + const int level = scheme->element_get_level (eclass, elements[0]); + + /* we set a maximum refinement level as forest user data */ + const int maxlevel = *(int *) t8_forest_get_user_data (forest); + if (level >= maxlevel) { + /* Do not refine after the maxlevel */ + return 0; + } + const t8_gloidx_t global_tree = t8_forest_global_tree_id (forest_from, which_tree); + if (global_tree == 0) { + return 1; + } + return 0; +} + +/** This test case loads a specific 2D .msh file and iterates over it to compute + * the leaf face neighbors. + * It originated from debugging session by Benedict Geihe using Trixi.jl where + * this setting resulted in errors. + * This, we implemented it as a test case in order to catch possible errors again in the future. + * The mesh consists of 2 quad trees that are connected via nontrivial orientation. + * + * - Load msh file with 2 quad trees + * - Adapt mesh such that the first tree is refined once uniformly. + * - Iterate over all mesh elements and compute the leaf face neighbors. + * + * __ __ _____ + * |__|__| | + * |__|__|_____| + * Tree_0 Tree_1 + * + * +*/ +class forest_face_neighbors_two_quad_mesh: public testing::TestWithParam { + protected: + void + SetUp () override + { + + /* Read our specific mesh file into a cmesh and build a forest. */ + const std::string meshfile_prefix = std::string (T8_TEST_DATA_DIR) + "/test_twosquares_twisted"; + const int partition_mesh = 0; + const int mesh_dim = 2; + const int main_proc = 0; + const int use_cad = 0; + t8_cmesh_t cmesh = t8_cmesh_from_msh_file (meshfile_prefix.c_str (), partition_mesh, sc_MPI_COMM_WORLD, mesh_dim, + main_proc, use_cad); + ASSERT_NE (cmesh, nullptr) << "Could not open mesh file."; + + /* Build a uniform forest on it. */ + const int scheme_id = GetParam (); + const t8_scheme *scheme = create_from_scheme_id (scheme_id); + const int level = 0; + const int do_ghost = 1; + t8_forest_t forest = t8_forest_new_uniform (cmesh, scheme, level, do_ghost, sc_MPI_COMM_WORLD); + + /* Build an adaptive forest */ + const int do_recursive_adapt = 0; + const int max_adapt_level = 2; + adaptive_forest + = t8_forest_new_adapt (forest, t8_test_adapt_first_tree, do_recursive_adapt, do_ghost, (void *) &max_adapt_level); + } + + void + TearDown () override + { + if (adaptive_forest != nullptr) { + t8_forest_unref (&adaptive_forest); + } + else if (cmesh != nullptr) { + t8_cmesh_unref (&cmesh); + } + } + + t8_cmesh_t cmesh = nullptr; + t8_forest_t adaptive_forest = nullptr; +}; + +/* Perform the actual test for the forest_face_neighbors_two_quad_mesh. + * Iterate over all leaves and ghosts and call the leaf face neighbor function. */ +TEST_P (forest_face_neighbors_two_quad_mesh, check_neighbors) +{ + // For debug purpoese, we write the forest to vtk + t8_forest_write_vtk (adaptive_forest, "lfn_test_twoquads"); + + // Iterate over all elements and compute their neighbors. + const t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (adaptive_forest); + const t8_locidx_t num_ghost_trees = t8_forest_get_num_ghost_trees (adaptive_forest); + t8_locidx_t ielement_index = 0; + for (t8_locidx_t itree = 0; itree < num_local_trees + num_ghost_trees; itree++) { + const t8_gloidx_t gtreeid = t8_forest_global_tree_id (adaptive_forest, itree); + const bool is_ghost = itree >= num_local_trees; + const t8_locidx_t ghost_tree_id = itree - num_local_trees; + /* Get the leaf element array */ + const t8_element_array_t *leaf_elements + = is_ghost ? t8_forest_ghost_get_tree_leaf_elements (adaptive_forest, ghost_tree_id) + : t8_forest_get_tree_leaf_element_array (adaptive_forest, itree); + const t8_eclass_t tree_class = t8_forest_get_tree_class (adaptive_forest, itree); + const t8_scheme *scheme = t8_forest_get_scheme (adaptive_forest); + const t8_locidx_t num_leaves = t8_element_array_get_count (leaf_elements); + int num_faces_with_2_neighbors = 0; + int num_faces_with_1_neighbor = 0; + for (t8_locidx_t ileaf = 0; ileaf < num_leaves; ++ileaf, ++ielement_index) { + // Iterate over each leaf element + const t8_element_t *element = t8_element_array_index_locidx (leaf_elements, ileaf); + const int num_faces = scheme->element_get_num_faces (tree_class, element); + for (int iface = 0; iface < num_faces; ++iface) { + // Iterate over all faces and compute the face neighbors + + // preparation + const t8_element_t **neighbor_leaves; + int *dual_faces; + int num_neighbors = 0; + t8_locidx_t *element_indices; + t8_eclass_t neigh_class; + t8_gloidx_t gneigh_tree; + int orientation; + + t8_debugf ("Compute face neighbor for tree %i (%s) element %i (index %i), at face %i.\n", itree, + is_ghost ? "ghost" : "local", ileaf, ielement_index, iface); + + // Actual computation of the face neighbors + t8_forest_leaf_face_neighbors_ext (adaptive_forest, itree, element, &neighbor_leaves, iface, &dual_faces, + &num_neighbors, &element_indices, &neigh_class, &gneigh_tree, &orientation); + + std::string buffer = ""; + for (int ineigh = 0; ineigh < num_neighbors; ++ineigh) { + buffer += std::to_string (element_indices[ineigh]) + " "; + } + t8_debugf ("Tree %i, Element %i, Face %i has %i neighbors:\t%s\n", itree, ileaf, iface, num_neighbors, + buffer.c_str ()); + + // We are now checking explicit face neighbor values. + // The number of neighbors must be >= 0 + EXPECT_GE (num_neighbors, 0); + // Count how many faces have 1 and 2 neighbors. + // We expect tree 1 to have 1 face with 2 neighbors + // We cannot state how many faces with 1 neighbor we count for tree 0, since we + // do not know how the elements are distributed across processes. + // But we know, that each element has at least one face with 1 neighbor, so the total count must be >0. + if (num_neighbors == 1) { + num_faces_with_1_neighbor++; + } + else if (num_neighbors == 2) { + num_faces_with_2_neighbors++; + } + if (gtreeid == 0) { + // Every face in global tree 0 has 0 or 1 face neighbors. + EXPECT_LE (num_neighbors, 1) << "Tree 0 faces must have 0 or 1 neighbor."; + } + else { + T8_ASSERT (gtreeid == 1); + // Tree 1 has only one element, it has either 0 neighbors, or 2 + if (num_neighbors > 0) { + EXPECT_EQ (num_neighbors, 2) << "Tree 2 faces must have 0 or 2 neighbors."; + } + } + + // clean-up original element neighbors + if (num_neighbors > 0) { + // Devnote: We are not using T8_TESTSUITE_FREE here, since the memory is allocated inside of t8_forest_leaf_face_neighbors_ext + // by t8code itself using T8_ALLOC. Thus, we need to call T8_FREE to free it. + T8_FREE (neighbor_leaves); + T8_FREE (element_indices); + T8_FREE (dual_faces); + } + } // End face loop + } // End leaf in tree loop + if (gtreeid == 0) { + // Tree 0 must have >0 faces with 1 neighbor. + EXPECT_GE (num_faces_with_1_neighbor, 1) << "Tree 1 must have 1 face with 2 neighbors."; + } + if (gtreeid == 1) { + // Tree 1 must have exactly 1 face with 2 neighbors. + EXPECT_EQ (num_faces_with_2_neighbors, 1) << "Tree 1 must have 1 face with 2 neighbors."; + } + } // End tree loop +} + +INSTANTIATE_TEST_SUITE_P (t8_gtest_face_neighbors, forest_face_neighbors_two_quad_mesh, AllSchemeCollections); diff --git a/test/testfiles/test_twosquares_twisted.msh b/test/testfiles/test_twosquares_twisted.msh new file mode 100644 index 0000000000..cb0920cf12 --- /dev/null +++ b/test/testfiles/test_twosquares_twisted.msh @@ -0,0 +1,84 @@ +$MeshFormat +4.1 0 8 +$EndMeshFormat +$Entities +6 7 2 0 +1 0 0 0 0 +2 1 0 0 0 +3 0 1 0 0 +4 1 1 0 0 +5 2 0 0 0 +6 2 1 0 0 +1 0 0 0 1 0 0 0 2 1 -2 +2 1 0 0 1 1 0 0 2 2 -4 +3 0 1 0 1 1 0 0 2 4 -3 +4 0 0 0 0 1 0 0 2 3 -1 +5 1 1 0 2 1 0 0 2 6 -4 +6 1 0 0 2 0 0 0 2 2 -5 +7 2 0 0 2 1 0 0 2 5 -6 +1 0 0 0 1 1 0 0 4 1 2 3 4 +2 1 0 0 2 1 0 0 4 5 -2 6 7 +$EndEntities +$Nodes +15 6 1 6 +0 1 0 1 +1 +0 0 0 +0 2 0 1 +2 +1 0 0 +0 3 0 1 +3 +0 1 0 +0 4 0 1 +4 +1 1 0 +0 5 0 1 +5 +2 0 0 +0 6 0 1 +6 +2 1 0 +1 1 0 0 +1 2 0 0 +1 3 0 0 +1 4 0 0 +1 5 0 0 +1 6 0 0 +1 7 0 0 +2 1 0 0 +2 2 0 0 +$EndNodes +$Elements +15 15 1 15 +0 1 15 1 +1 1 +0 2 15 1 +2 2 +0 3 15 1 +3 3 +0 4 15 1 +4 4 +0 5 15 1 +5 5 +0 6 15 1 +6 6 +1 1 1 1 +7 1 2 +1 2 1 1 +8 2 4 +1 3 1 1 +9 4 3 +1 4 1 1 +10 3 1 +1 5 1 1 +11 6 4 +1 6 1 1 +12 2 5 +1 7 1 1 +13 5 6 +2 1 3 1 +14 1 2 4 3 +2 2 3 1 +15 6 4 2 5 +$EndElements diff --git a/tutorials/general/t8_step6_stencil.cxx b/tutorials/general/t8_step6_stencil.cxx index ab5908fa2c..4f55bdfa18 100644 --- a/tutorials/general/t8_step6_stencil.cxx +++ b/tutorials/general/t8_step6_stencil.cxx @@ -206,15 +206,15 @@ t8_step6_compute_stencil (t8_forest_t forest, struct data_per_element *element_d /* Loop over all faces of an element. */ int num_faces = scheme->element_get_num_faces (tree_class, element); for (int iface = 0; iface < num_faces; iface++) { - int num_neighbors; /**< Number of neighbors for each face */ - int *dual_faces; /**< The face indices of the neighbor elements */ - t8_locidx_t *neighids; /**< Indices of the neighbor elements */ - t8_element_t **neighbors; /*< Neighboring elements. */ - t8_eclass_t neigh_class; /*< Neighboring elements tree class. */ + int num_neighbors; /**< Number of neighbors for each face */ + int *dual_faces; /**< The face indices of the neighbor elements */ + t8_locidx_t *neighids; /**< Indices of the neighbor elements */ + const t8_element_t **neighbors; /*< Neighboring elements. */ + t8_eclass_t neigh_class; /*< Neighboring elements tree class. */ /* Collect all neighbors at the current face. */ t8_forest_leaf_face_neighbors (forest, itree, element, &neighbors, iface, &dual_faces, &num_neighbors, - &neighids, &neigh_class, 1); + &neighids, &neigh_class); /* Retrieve the `height` of the face neighbor. Account for two neighbors in case of a non-conforming interface by computing the average. */ @@ -248,7 +248,6 @@ t8_step6_compute_stencil (t8_forest_t forest, struct data_per_element *element_d if (num_neighbors > 0) { /* Free allocated memory. */ - scheme->element_destroy (tree_class, num_neighbors, neighbors); T8_FREE (neighbors); T8_FREE (dual_faces); T8_FREE (neighids);