diff --git a/example/advect/t8_advection.cxx b/example/advect/t8_advection.cxx index 2ba506c9dd..866fb918ad 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 6b987f871c..8b3180c560 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_tree_get_leaf_elements (forest, itree) + : t8_forest_ghost_get_tree_leaf_elements (forest, ghost_tree_id); + // 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 5cf3ec1d8e..398d793331 100644 --- a/mesh_handle/element.hxx +++ b/mesh_handle/element.hxx @@ -306,14 +306,14 @@ class element: public TCompetences>... { } } std::vector> neighbor_elements; - 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, get_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); } @@ -329,7 +329,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_data/t8_containers.cxx b/src/t8_data/t8_containers.cxx index e025c9ed6b..fe20b2864c 100644 --- a/src/t8_data/t8_containers.cxx +++ b/src/t8_data/t8_containers.cxx @@ -27,6 +27,7 @@ #include #include #include +#include T8_EXTERN_C_BEGIN (); @@ -326,6 +327,39 @@ t8_element_array_get_array_mutable (t8_element_array_t *element_array) return &element_array->array; } +t8_locidx_t +t8_element_array_find (const t8_element_array_t *element_array, const t8_element_t *element) +{ + /* In order to find the element, we need to compute its linear id. + * To do so, we need the scheme and the level of the element. */ + const t8_scheme *scheme = t8_element_array_get_scheme (element_array); + const t8_eclass_t tree_class = t8_element_array_get_tree_class (element_array); + const int element_level = scheme->element_get_level (tree_class, element); + /* Compute the linear id. */ + const t8_linearidx_t element_id = scheme->element_get_linear_id (tree_class, element, element_level); + /* Search for the element. + * The search returns the largest index i, + * such that the element at position i has a smaller id than the given one. + * If no such i exists, it returns -1. */ + const t8_locidx_t search_result = t8_forest_bin_search_lower (element_array, element_id, element_level); + if (search_result < 0) { + // The element was not found, we return -1. */ + return -1; + } + /* An element was found but it may not be the candidate element. + * To identify whether the element was found, we compare these two. */ + const t8_element_t *check_element = t8_element_array_index_locidx (element_array, search_result); + T8_ASSERT (check_element != NULL); + if (scheme->element_is_equal (tree_class, element, check_element)) { + // The element was found at position search_result. We return it. + return search_result; + } + else { + // The element was not found, we return -1. */ + return -1; + } +} + void t8_element_array_reset (t8_element_array_t *element_array) { diff --git a/src/t8_data/t8_containers.h b/src/t8_data/t8_containers.h index efa2ff6d6b..f79e79b8f6 100644 --- a/src/t8_data/t8_containers.h +++ b/src/t8_data/t8_containers.h @@ -265,6 +265,16 @@ t8_element_array_get_array (const t8_element_array_t *element_array); sc_array_t * t8_element_array_get_array_mutable (t8_element_array_t *element_array); +/** Search for an element in an array. + * \param [in] element_array Array structure. + * \param [in] element Element to be found in \a element_array. + * The element must have been created with the scheme used in \a element_array. + * \return If \a element was found in \a element_array then the position in the array is returned. + * If the element is not found, -1 is returned. +*/ +t8_locidx_t +t8_element_array_find (const t8_element_array_t *element_array, const t8_element_t *element); + /** Sets the array count to zero and frees all elements. * \param [in,out] element_array Array structure to be reset. * \note Calling t8_element_array_init, then any array operations, diff --git a/src/t8_forest/t8_forest.cxx b/src/t8_forest/t8_forest.cxx index c02b2e52ce..152d0fa67f 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).*/ 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: If the forest has no ghosts, then skip the ghosts + parts. In that case, process boundary elements will have 0 neighbors. +*/ t8_gloidx_t t8_forest_element_face_neighbor (t8_forest_t forest, t8_locidx_t ltreeid, const t8_element_t *elem, t8_element_t *neigh, - 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 lctree_id = t8_forest_ltreeid_to_cmesh_ltreeid (forest, ltreeid); + T8_ASSERT (lctree_id >= 0); +#if T8_ENABLE_DEBUG + const bool cmesh_tree_is_local = t8_cmesh_treeid_is_local_tree (cmesh, lctree_id); + T8_ASSERT (cmesh_tree_is_local || t8_cmesh_treeid_is_ghost (cmesh, lctree_id)); +#endif + + const t8_locidx_t neighbor_ctreeid + = t8_cmesh_get_face_neighbor (cmesh, lctree_id, tree_face, &tree_neigh_face, &orientation); + + if (neighbor_ctreeid < 0) { /* This face is a domain boundary. We do not need to continue */ return -1; } + /* Get the eclass for the boundary */ const t8_eclass_t boundary_class = (t8_eclass_t) t8_eclass_face_types[eclass][tree_face]; /* Allocate the face element */ scheme->element_new (boundary_class, 1, &face_element); /* Compute the face element. */ scheme->element_get_boundary_face (eclass, elem, face, face_element); - /* 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,473 @@ 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 : ltreeid - t8_forest_get_num_local_trees (forest); + T8_ASSERT (t8_forest_element_is_leaf_or_ghost (forest, element, adjusted_tree_id, is_ghost_tree)); + + struct t8_lfn_user_data *lfn_data = reinterpret_cast (user_data); + // face is the face of the considered leaf neighbor element and thus the + // corresponding dual face + t8_debugf ("Adding new face neighbor (leaf index %i) with dual face %i.\n", tree_leaf_index, face); + lfn_data->dual_faces.push_back (face); + // Compute the index of the element + const t8_locidx_t num_local_elements = t8_forest_get_local_num_leaf_elements (forest); + const t8_locidx_t tree_offset + = !is_ghost_tree ? t8_forest_get_tree_element_offset (forest, ltreeid) + : t8_forest_ghost_get_tree_element_offset (forest, adjusted_tree_id) + num_local_elements; + const t8_locidx_t element_index = tree_offset + tree_leaf_index; + lfn_data->element_indices.push_back (element_index); + // Add the pointer to the current element + const t8_element_t *&pnew_element = lfn_data->neighbors.emplace_back (); + if (!is_ghost_tree) { + pnew_element = t8_forest_get_leaf_element_in_tree (forest, ltreeid, tree_leaf_index); + } + else { + pnew_element = t8_forest_ghost_get_leaf_element (forest, adjusted_tree_id, tree_leaf_index); + } + return 1; +} + +/* + Set the proper return values of the leaf face neighbor computation + in case that no neighbors are found. + */ +static void +t8_forest_leaf_face_neighbors_set_no_neighbor_return_value (const t8_element_t **pneighbor_leaves[], int *dual_faces[], + int *num_neighbors, t8_locidx_t **pelement_indices, + t8_gloidx_t *gneigh_tree) +{ + *dual_faces = NULL; + *num_neighbors = 0; + *pelement_indices = NULL; + if (pneighbor_leaves) { + *pneighbor_leaves = NULL; + } + if (gneigh_tree != NULL) { + *gneigh_tree = -1; + } +} + +/* TODO: If the forest has no ghosts, then skip the ghosts + parts. In that case, process boundary elements will have 0 neighbors. +*/ +void +t8_forest_leaf_face_neighbors_ext (const t8_forest_t forest, const t8_locidx_t ltreeid, + const t8_element_t *leaf_or_ghost, const t8_element_t **pneighbor_leaves[], + const int face, int *dual_faces[], int *num_neighbors, + t8_locidx_t **pelement_indices, t8_eclass_t *pneigh_eclass, t8_gloidx_t *gneigh_tree, + int *orientation) +{ + /* We compute all face neighbor leaf elements of E via the following strategy: + * - Compute the same level face neighbor N + * - The neighbor tree could be a local tree or ghost (or both), + * for each variant get the leaf array of the neighbor tree and search in it: + * - Search for the first leaf element L overlapping with N. + * If it exists, it is either an ancestor or descendant of N (or N itself which is included in both definitions). + * If it does not exist, there are not leaf face neighbors in this tree. + * - If L is an ancestor of N (i.e. Level(L) <= Level(N)) then L is the only face neighbor. + * - Otherwise (Level(L) > Level (N)) we use a recursive face search across N's neighbor face, + * adding all leaf elements on the face to the face neighbors. + * This search will require L (more precise its position in the tree leaf array) as input. + **/ T8_ASSERT (t8_forest_is_committed (forest)); - T8_ASSERT (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 +2103,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 +2113,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); @@ -1957,36 +2140,36 @@ t8_forest_tree_is_local (const t8_forest_t forest, const t8_locidx_t local_tree) int t8_forest_element_is_leaf (const t8_forest_t forest, const t8_element_t *element, const t8_locidx_t local_tree) { - T8_ASSERT (t8_forest_is_committed (forest)); + bool check_ghost = false; T8_ASSERT (t8_forest_tree_is_local (forest, local_tree)); + return t8_forest_element_is_leaf_or_ghost (forest, element, local_tree, check_ghost); +} + +int +t8_forest_element_is_leaf_or_ghost (const t8_forest_t forest, const t8_element_t *element, const t8_locidx_t local_tree, + const int check_ghost) +{ + T8_ASSERT (t8_forest_is_committed (forest)); +#if T8_ENABLE_DEBUG + if (!check_ghost) { + T8_ASSERT (t8_forest_tree_is_local (forest, local_tree)); + } + else { + T8_ASSERT (0 <= local_tree && local_tree < t8_forest_get_num_ghost_trees (forest)); + } +#endif /* We get the array of the tree's elements and then search in the array of elements for our * element candidate. */ /* Get the array */ - const t8_element_array_t *elements = t8_forest_get_tree_leaf_element_array (forest, local_tree); + const t8_element_array_t *elements = !check_ghost ? t8_forest_get_tree_leaf_element_array (forest, local_tree) + : t8_forest_ghost_get_tree_leaf_elements (forest, local_tree); + T8_ASSERT (elements != NULL); - /* In order to find the element, we need to compute its linear id. - * To do so, we need the scheme and the level of the element. */ - const t8_scheme *scheme = t8_element_array_get_scheme (elements); - const t8_eclass_t tree_class = t8_element_array_get_tree_class (elements); - const int element_level = scheme->element_get_level (tree_class, element); - /* Compute the linear id. */ - const t8_linearidx_t element_id = scheme->element_get_linear_id (tree_class, element, element_level); - /* Search for the element. - * The search returns the largest index i, - * such that the element at position i has a smaller id than the given one. - * If no such i exists, it returns -1. */ - const t8_locidx_t search_result = t8_forest_bin_search_lower (elements, element_id, element_level); - if (search_result < 0) { - /* The element was not found. */ - return 0; - } - /* An element was found but it may not be the candidate element. - * To identify whether the element was found, we compare these two. */ - const t8_element_t *check_element = t8_element_array_index_locidx (elements, search_result); - T8_ASSERT (check_element != NULL); - return (scheme->element_is_equal (tree_class, element, check_element)); + // Search for the element in the array, return true if it was found, + // false if not. + return t8_element_array_find (elements, element) >= 0; } /* Check if an element is owned by a specific rank */ @@ -2555,7 +2738,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 +2771,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 a6d171d187..7c985b7437 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 @@ -550,6 +550,25 @@ t8_forest_get_coarse_tree (t8_forest_t forest, t8_locidx_t ltreeid); int t8_forest_element_is_leaf (const t8_forest_t forest, const t8_element_t *element, const t8_locidx_t local_tree); +/** + * Query whether a given element or a ghost is a leaf of a local or ghost tree in a forest. + * + * \param [in] forest The forest. + * \param [in] element An element of a local tree in \a forest. + * \param [in] local_tree A local tree id of \a forest or a ghost tree id + * \param [in] check_ghost If true \a element is interpreted as a ghost element and + * \a local_tree as the id of a ghost tree (0 <= \a local_tree < num_ghost_trees). + * If false \a element is interpreted as an element and \a local_tree as + * the id of a local tree (0 <= \a local_tree < num_local_trees). + * \return True (non-zero) if and only if \a element is a leaf (or ghost) in \a local_tree of \a forest. + * \note \a forest must be committed before calling this function. + * \ref t8_forest_element_is_leaf + * \ref t8_forest_element_is_ghost + */ +int +t8_forest_element_is_leaf_or_ghost (const t8_forest_t forest, const t8_element_t *element, const t8_locidx_t local_tree, + const int check_ghost); + /** Compute the leaf face orientation at given face in a forest. * \param [in] forest The forest. Must have a valid ghost layer. * \param [in] ltreeid A local tree id. @@ -565,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. @@ -580,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. (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); @@ -598,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 @@ -617,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); @@ -640,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. @@ -836,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. @@ -863,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 5e600bf746..27fcb22874 100644 --- a/src/t8_forest/t8_forest_ghost.cxx +++ b/src/t8_forest/t8_forest_ghost.cxx @@ -337,6 +337,16 @@ t8_forest_ghost_num_trees (const t8_forest_t forest) return forest->ghosts->ghost_trees->elem_count; } +#if T8_ENABLE_DEBUG +static bool +t8_forest_tree_is_ghost (const t8_forest_t forest, const t8_locidx_t lghost_tree) +{ + T8_ASSERT (t8_forest_is_committed (forest)); + + return 0 <= lghost_tree && lghost_tree < t8_forest_get_num_ghost_trees (forest); +} +#endif + /** * Given an index into the ghost_trees array return the ghost tree * @@ -355,7 +365,7 @@ t8_forest_ghost_get_tree (const t8_forest_t forest, const t8_locidx_t lghost_tre ghost = forest->ghosts; T8_ASSERT (ghost != NULL); T8_ASSERT (ghost->ghost_trees != NULL); - T8_ASSERT (0 <= lghost_tree && lghost_tree < t8_forest_ghost_num_trees (forest)); + T8_ASSERT (t8_forest_tree_is_ghost (forest, lghost_tree)); ghost_tree = (t8_ghost_tree_t *) t8_sc_array_index_locidx (ghost->ghost_trees, lghost_tree); return ghost_tree; @@ -446,6 +456,14 @@ t8_forest_ghost_get_leaf_element (t8_forest_t forest, t8_locidx_t lghost_tree, t return t8_element_array_index_locidx_mutable (&ghost_tree->elements, lelement); } +int +t8_forest_element_is_ghost (const t8_forest_t forest, const t8_element_t *element, const t8_locidx_t lghost_tree) +{ + bool check_ghost = true; + T8_ASSERT (t8_forest_tree_is_ghost (forest, lghost_tree)); + return t8_forest_element_is_leaf_or_ghost (forest, element, lghost_tree, check_ghost); +} + /** Initialize a t8_ghost_remote_tree_t. * * \param[in] forest The forest. @@ -846,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 c00c39a1af..cf36bbf51e 100644 --- a/src/t8_forest/t8_forest_ghost.h +++ b/src/t8_forest/t8_forest_ghost.h @@ -86,7 +86,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. */ @@ -118,7 +118,7 @@ t8_forest_ghost_get_tree_class (const t8_forest_t forest, const t8_locidx_t lgho /** Given a local ghost tree compute the global tree id of it. * \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 The global id of the local ghost tree \a lghost_tree. * \a forest must be committed before calling this function. * \see https://github.com/DLR-AMR/t8code/wiki/Tree-indexing for more details about tree indexing. @@ -137,6 +137,18 @@ t8_forest_ghost_get_global_treeid (const t8_forest_t forest, const t8_locidx_t l t8_element_t * t8_forest_ghost_get_leaf_element (t8_forest_t forest, t8_locidx_t lghost_tree, t8_locidx_t lelement); +/** + * Query whether a given element is a ghost of a certrain tree in a forest. + * + * \param [in] forest The forest. + * \param [in] element An element of a ghost tree in \a forest. + * \param [in] lghost_tree A local ghost tree id of \a forest. (0 <= \a lghost_tree < num_ghost_trees) + * \return True (non-zero) if and only if \a element is a ghost in \a lghost_tree of \a forest. + * \note \a forest must be committed before calling this function. + */ +int +t8_forest_element_is_ghost (const t8_forest_t forest, const t8_element_t *element, const t8_locidx_t lghost_tree); + /** Return the array of remote ranks. * \param [in] forest A forest with constructed ghost layer. * \param [in,out] num_remotes On output the number of remote ranks is stored here. 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 deb30c3c60..f10e3698ea 100644 --- a/src/t8_forest/t8_forest_iterate.h +++ b/src/t8_forest/t8_forest_iterate.h @@ -40,14 +40,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 @@ -152,6 +155,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. @@ -166,14 +173,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.cxx b/src/t8_forest/t8_forest_private.cxx index ce3e91879d..49d3367eda 100644 --- a/src/t8_forest/t8_forest_private.cxx +++ b/src/t8_forest/t8_forest_private.cxx @@ -57,6 +57,14 @@ t8_forest_get_tree_leaf_element_array_mutable (const t8_forest_t forest, t8_loci return (t8_element_array_t *) t8_forest_get_tree_leaf_element_array (forest, ltreeid); } +/* TODO: does the search fail when element_level is smaller then levels in the array? + For example entering the search with the root element or a level 1 element + and the array contains much finer elements. + Will it still return the largest index, or just any index? + */ +/* TODO: This may be implementable with std::partition_point, which would yield an easier implementation. + Need to check. + */ /** \brief Search for a linear element id (at level element_level) in a sorted array of * elements. If the element does not exist, return the largest index i * such that the element at position i has a smaller id than the given one. @@ -93,4 +101,123 @@ t8_forest_bin_search_lower (const t8_element_array_t *elements, const t8_lineari return elem_iter.get_current_index () - 1; } +/** \brief Search for a linear element id (at level element_level) in a sorted array of + * elements. If the element does not exist, return the smallest index i + * such that the element at position i has a larger id than the given one. + * If no such i exists, return -1. + */ +t8_locidx_t +t8_forest_bin_search_upper (const t8_element_array_t *elements, const t8_linearidx_t element_id, + const int element_level) +{ + const t8_scheme *scheme = t8_element_array_get_scheme (elements); + const t8_eclass_t tree_class = t8_element_array_get_tree_class (elements); + /* At first, we check whether any element has smaller id than the + * given one. */ + const t8_locidx_t num_elements = t8_element_array_get_count (elements); + if (num_elements == 0) { + /* This array is empty. */ + return -1; + } + const t8_element_t *query = t8_element_array_index_int (elements, num_elements - 1); + const t8_linearidx_t query_id = scheme->element_get_linear_id (tree_class, query, element_level); + if (query_id < element_id) { + /* No element has id larger than the given one. */ + return -1; + } + + /* We search for the first element E in the array, where element_id > ID(E) is false. + Thus, E is the first element with ID(E) >= element_id . */ + auto elem_iter + = std::lower_bound (t8_element_array_begin (elements), t8_element_array_end (elements), element_id, + [&element_level, &scheme, &tree_class] (const t8_element_array_iterator::value_type &elem_ptr, + const t8_linearidx_t element_id_) { + return (element_id_ > scheme->element_get_linear_id (tree_class, elem_ptr, element_level)); + }); + + /* In case we do not find an element that is greater than the given element_id, the binary search returns + * the end-iterator of the element array. */ + if (elem_iter == t8_element_array_end (elements)) { + // No element was found. + return -1; + } + else { + return elem_iter.get_current_index (); + } +} + +/** \brief Search for a linear element id (at level element_level) in a sorted array of + * elements. If the element does not exist, return the first index i such that + * the element at position i is an ancestor or descendant of the element corresponding to the element id. + * If no such i exists, return -1. + */ +t8_locidx_t +t8_forest_bin_search_first_descendant_ancestor (const t8_element_array_t *elements, const t8_element_t *element, + const t8_element_t **element_found) +{ + /* This search works as follows: + + Let E denote the element with element_id at level L. + If an ancestor or descendant of E exists in the array then they are either: + A: The search result of t8_forest_bin_search_lower + B: The search result of t8_forest_bin_search_upper + + Let ID(element,level) denote the linear id of an element at a given level. + + Case A: There is an element F in the array that is E itself or an ancestor of E (i.e. level(F) <= level(E)). + In that case + ID(E,L) >= ID(F,L) and there can be no element with id in between (since it would also be an ancestor of E). + Then F will be the search result of t8_forest_bin_search_lower + Case B: There is an element F in the array that is a descendant of E and it has the smallest index in the array of all descendants. + Then + ID(E,L) = ID(F,L) + and also + ID(E,L) = ID(D,L) for all other descendants of E. + But since F is the first it will be the search result of t8_forest_bin_search_upper. + Case C: There is no descendant or ancestor of E in the array. In both cases t8_forest_bin_search_lower and + t8_forest_bin_search_upper may find elements but the results will not be ancestors/descendants of E. + + From this, we determine the following algorithm: + + 1. Query t8_forest_bin_search_lower with N. + 2. If no element was found, or the resulting element is not an ancestor of N. + 3. Query t8_forest_bin_search_upper with N. + 3. If an element was found and it is a descendant of N, we found our element. + 4. If not, no element was found. + */ + + /* Compute the element's level and linear id. In order to do so, + * we first need the scheme and eclass. */ + const t8_scheme *scheme = t8_element_array_get_scheme (elements); + const t8_eclass eclass = t8_element_array_get_tree_class (elements); + const int element_level = scheme->element_get_level (eclass, element); + const t8_linearidx_t element_id = scheme->element_get_linear_id (eclass, element, element_level); + + const t8_locidx_t search_pos_lower = t8_forest_bin_search_lower (elements, element_id, element_level); + + /* Get the element at the current position. */ + if (search_pos_lower >= 0) { + *element_found = t8_element_array_index_locidx (elements, search_pos_lower); + const bool is_ancestor = scheme->element_is_ancestor (eclass, *element_found, element); + if (is_ancestor) { + /* The element at this position is an ancestor or descendant. */ + return search_pos_lower; + } + } + /* t8_forest_bin_search_lower did not return a result or an ancestor. */ + + const t8_locidx_t search_pos_upper = t8_forest_bin_search_upper (elements, element_id, element_level); + if (search_pos_upper >= 0) { + *element_found = t8_element_array_index_locidx (elements, search_pos_upper); + const bool is_descendant = scheme->element_is_ancestor (eclass, element, *element_found); + if (is_descendant) { + /* The element at this position is an ancestor or descendant. */ + return search_pos_upper; + } + } + // No ancestor or descendant was found + *element_found = nullptr; + return -1; +} + T8_EXTERN_C_END (); diff --git a/src/t8_forest/t8_forest_private.h b/src/t8_forest/t8_forest_private.h index 84c8f485a4..d872d8023b 100644 --- a/src/t8_forest/t8_forest_private.h +++ b/src/t8_forest/t8_forest_private.h @@ -217,6 +217,35 @@ t8_locidx_t t8_forest_bin_search_lower (const t8_element_array_t *elements, const t8_linearidx_t element_id, const int element_level); +/** \brief Search for a linear element id (at level element_level) in a sorted array of + * elements. If the element does not exist, return the smallest index i + * such that the element at position i has a larger id than the given one. + * If no such i exists, return -1. + * \param [in] elements An array of elements. Must be sorted according to linear id at maximum level. + * Must correspond to a valid refinement (i.e. contain no duplicate elements or elements and their descendants). + * \param [in] element_id The linear id of the element to search for. + * \param [in] element_level The level of the element to search for. Thus, the level at which \a element_id was computed. + * \return The smallest index \a i of an element with linear_id larger than or equal to \a element_id in \a elements if it exists. + * -1 if no such element was found in \a elements. + */ +t8_locidx_t +t8_forest_bin_search_upper (const t8_element_array_t *elements, const t8_linearidx_t element_id, + const int element_level); + +/** \brief Search for the first descendant or ancestor of an element in a sorted array of elements. + * \param [in] elements An array of elements. Must be sorted according to linear id at maximum level. + * Must correspond to a valid refinement (i.e. contain no duplicate elements or elements and their descendants). + * \param [in] element The element to search for. + * \param [in] element_found On return either a descendant or ancestor of \a element in \a elements if it exists. NULL if no + * such element exists in \a elements. + * \return The smallest index \a i such that elements[i] (= \a element_found) is an ancestor or a descendant of \a element. + * -1 if no such element was found in \a elements. + * \note \a element is ancestor and descendant of itself, so if \a element is contained in \a elements then it will be found by this function. + */ +t8_locidx_t +t8_forest_bin_search_first_descendant_ancestor (const t8_element_array_t *elements, const t8_element_t *element, + const t8_element_t **element_found); + /** Find the owner process of a given element, deprecated version. * Use t8_forest_element_find_owner instead. * \param [in] forest The forest. @@ -337,7 +366,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/src/t8_schemes/t8_default/t8_default_common/t8_default_common.hxx b/src/t8_schemes/t8_default/t8_default_common/t8_default_common.hxx index f1f8146ff2..169430dd49 100644 --- a/src/t8_schemes/t8_default/t8_default_common/t8_default_common.hxx +++ b/src/t8_schemes/t8_default/t8_default_common/t8_default_common.hxx @@ -191,6 +191,40 @@ struct t8_default_scheme_common: public t8_scheme_helpersunderlying ().element_is_valid (element_A)); + T8_ASSERT (this->underlying ().element_is_valid (element_B)); + + const int level_A = this->underlying ().element_get_level (element_A); + const int level_B = this->underlying ().element_get_level (element_B); + + if (level_A > level_B) { + /* element A is finer than element B and thus cannot be + * an ancestor of B. */ + return false; + } + + const t8_linearidx_t id_A = this->underlying ().element_get_linear_id (element_A, level_A); + const t8_linearidx_t id_B = this->underlying ().element_get_linear_id (element_B, level_A); + + // If both elements have the same linear ID at level_A then A is an ancestor of B. + return id_A == id_B; + } + /** Allocate space for a bunch of elements. * \param [in] length The number of elements to allocate. * \param [out] elem The elements to allocate. diff --git a/src/t8_schemes/t8_scheme.hxx b/src/t8_schemes/t8_scheme.hxx index fd398584a7..4124ca9af5 100644 --- a/src/t8_schemes/t8_scheme.hxx +++ b/src/t8_schemes/t8_scheme.hxx @@ -509,6 +509,21 @@ struct t8_scheme eclass_schemes[tree_class]); }; + /** Query whether element A is an ancestor of the element B. + * An element A is ancestor of an element B if A == B or if B can + * be obtained from A via successive refinement. + * \param [in] tree_class The eclass of the current tree. + * \param [in] element_A An element of class \a eclass in scheme \a scheme. + * \param [in] element_B An element of class \a eclass in scheme \a scheme. + * \return True if and only if \a element_A is an ancestor of \a element_B. + */ + bool + element_is_ancestor (const t8_eclass_t tree_class, const t8_element_t *element_A, const t8_element_t *element_B) const + { + return std::visit ([&] (auto &&scheme) { return scheme.element_is_ancestor (element_A, element_B); }, + eclass_schemes[tree_class]); + } + /** Query whether a given set of elements is a family or not. * \param [in] tree_class The eclass of the current tree. * \param [in] fam An array of as many elements as an element of class @@ -537,8 +552,9 @@ struct t8_scheme element_get_nca (const t8_eclass_t tree_class, const t8_element_t *elem1, const t8_element_t *elem2, t8_element_t *const nca) const { - return std::visit ([&] (auto &&scheme) { return scheme.element_get_nca (elem1, elem2, nca); }, - eclass_schemes[tree_class]); + std::visit ([&] (auto &&scheme) { return scheme.element_get_nca (elem1, elem2, nca); }, eclass_schemes[tree_class]); + T8_ASSERT (element_is_ancestor (tree_class, nca, elem1)); + T8_ASSERT (element_is_ancestor (tree_class, nca, elem2)); }; /** Compute the shape of the face of an element. diff --git a/src/t8_schemes/t8_standalone/t8_standalone_implementation.hxx b/src/t8_schemes/t8_standalone/t8_standalone_implementation.hxx index 5338e72625..9817224acd 100644 --- a/src/t8_schemes/t8_standalone/t8_standalone_implementation.hxx +++ b/src/t8_schemes/t8_standalone/t8_standalone_implementation.hxx @@ -608,6 +608,54 @@ struct t8_standalone_scheme: public t8_scheme_helpers level(B) then A cannot be an ancestor. + - Otherwise compute the ancestor of B at level(A) + - Compare the computed ancestor with A. + */ + T8_ASSERT (element_is_valid (element_A)); + T8_ASSERT (element_is_valid (element_B)); + + const t8_standalone_element *el_B = (const t8_standalone_element *) element_B; + + const int level_A = element_get_level (element_A); + const int level_B = element_get_level (element_B); + + if (level_A > level_B) { + // A is finer than B and thus cannot be an ancestor. + return false; + } + + // Compute the ancestor of B at level_A and compare it with A + t8_element_t *ancestor; + element_new (1, &ancestor); + + t8_standalone_element *ancestor_casted = (t8_standalone_element *) ancestor; + + element_get_ancestor (el_B, level_A, ancestor_casted); + + const bool is_ancestor = element_is_equal (ancestor, element_A); + + element_destroy (1, &ancestor); + + // Return true if A == ancestor + // Return false if A != ancestor + return is_ancestor; + } + /** Compute the nearest common ancestor of two elements. That is, * the element with highest level that still has both given elements as * descendants. diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6d28bbadb8..701f481095 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -149,7 +149,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_element_is_leaf_serial SOURCES t8_forest/t8_gtest_element_is_leaf.cxx ) +add_t8_cpp_test( NAME t8_gtest_element_is_leaf_serial SOURCES t8_forest/t8_gtest_element_is_leaf.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 ) @@ -184,6 +184,7 @@ add_t8_cpp_test( NAME t8_gtest_pyra_connectivity_serial SOURCES t8_schemes/t add_t8_cpp_test( NAME t8_gtest_face_neigh_serial SOURCES t8_schemes/t8_gtest_face_neigh.cxx ) add_t8_cpp_test( NAME t8_gtest_get_linear_id_serial SOURCES t8_schemes/t8_gtest_get_linear_id.cxx ) add_t8_cpp_test( NAME t8_gtest_ancestor_serial SOURCES t8_schemes/t8_gtest_ancestor.cxx ) +add_t8_cpp_test( NAME t8_gtest_is_ancestor_serial SOURCES t8_schemes/t8_gtest_is_ancestor.cxx ) add_t8_cpp_test( NAME t8_gtest_ancestor_id_serial SOURCES t8_schemes/t8_gtest_ancestor_id.cxx ) add_t8_cpp_test( NAME t8_gtest_element_count_leaves_serial SOURCES t8_schemes/t8_gtest_element_count_leaves.cxx ) add_t8_cpp_test( NAME t8_gtest_element_ref_coords_serial SOURCES t8_schemes/t8_gtest_element_ref_coords.cxx ) diff --git a/test/mesh_handle/t8_gtest_ghost.cxx b/test/mesh_handle/t8_gtest_ghost.cxx index bc5e600b9d..a31015b0f0 100644 --- a/test/mesh_handle/t8_gtest_ghost.cxx +++ b/test/mesh_handle/t8_gtest_ghost.cxx @@ -116,14 +116,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); @@ -139,7 +138,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 f19af2371c..11d10c6890 100644 --- a/test/t8_forest/t8_gtest_element_is_leaf.cxx +++ b/test/t8_forest/t8_gtest_element_is_leaf.cxx @@ -25,9 +25,11 @@ #include #include #include +#include #include #include "test/t8_cmesh_generator/t8_cmesh_example_sets.hxx" #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. @@ -42,34 +44,8 @@ #else #define T8_IS_LEAF_MAX_LVL 4 #endif -/* Adapt a forest such that always the first child of a - * family is refined and no other elements. This results in a highly - * imbalanced forest. */ -static int -t8_test_adapt_first_child (t8_forest_t forest, [[maybe_unused]] t8_forest_t forest_from, - [[maybe_unused]] t8_locidx_t which_tree, const t8_eclass_t tree_class, - [[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 (tree_class, elements[0]))); - - const int level = scheme->element_get_level (tree_class, elements[0]); - - /* we set a maximum refinement level as forest user data */ - int maxlevel = *(int *) t8_forest_get_user_data (forest); - if (level >= maxlevel) { - /* Do not refine after the maxlevel */ - return 0; - } - const int child_id = scheme->element_get_child_id (tree_class, elements[0]); - if (child_id == 1) { - return 1; - } - return 0; -} -struct element_is_leaf: public testing::TestWithParam, int>> +struct element_is_leaf_or_ghost: public testing::TestWithParam, int>> { protected: void @@ -82,9 +58,8 @@ struct element_is_leaf: public testing::TestWithParam (GetParam ()); t8_cmesh_t cmesh = t8_cmesh_new_from_class (tree_class, sc_MPI_COMM_WORLD); - forest = t8_forest_new_uniform (cmesh, scheme, level, 0, sc_MPI_COMM_WORLD); + forest = t8_forest_new_uniform (cmesh, scheme, level, 1, sc_MPI_COMM_WORLD); t8_forest_ref (forest); - //const int maxlevel = t8_forest_get_maxlevel (forest); int maxlevel = 7; const int recursive_adapt = 1; forest_adapt = t8_forest_new_adapt (forest, t8_test_adapt_first_child, recursive_adapt, 0, &maxlevel); @@ -106,7 +81,7 @@ struct element_is_leaf: public testing::TestWithParam +struct element_is_leaf_or_ghost_hybrid: public testing::TestWithParam { protected: void @@ -153,48 +128,97 @@ t8_test_element_is_leaf_for_forest (t8_forest_t forest) t8_element_t *not_leaf; scheme->element_new (tree_class, 1, ¬_leaf); /* Iterate over all the tree's leaf elements, check whether the leaf - * is correctly identified by t8_forest_element_is_leaf, + * is correctly identified by t8_forest_element_is_leaf and t8_forest_element_is_leaf_or_ghost, * build its parent and its first child (if they exist), and verify - * that t8_forest_element_is_leaf returns false. */ + * that t8_forest_element_is_leaf and t8_forest_element_is_leaf_or_ghost returns false. */ for (t8_locidx_t ielement = 0; ielement < num_elements_in_tree; ++ielement) { const t8_element_t *leaf_element = t8_forest_get_leaf_element_in_tree (forest, itree, ielement); EXPECT_TRUE (t8_forest_element_is_leaf (forest, leaf_element, itree)); + EXPECT_TRUE (t8_forest_element_is_leaf_or_ghost (forest, leaf_element, itree, 0)); /* Compute parent and first child of element and check that they are not in the tree */ const int element_level = scheme->element_get_level (tree_class, leaf_element); if (element_level > 0) { scheme->element_get_parent (tree_class, leaf_element, not_leaf); EXPECT_FALSE (t8_forest_element_is_leaf (forest, not_leaf, itree)); + EXPECT_FALSE (t8_forest_element_is_leaf_or_ghost (forest, not_leaf, itree, 0)); } if (element_level < scheme->get_maxlevel (tree_class)) { scheme->element_get_child (tree_class, leaf_element, 0, not_leaf); EXPECT_FALSE (t8_forest_element_is_leaf (forest, not_leaf, itree)); + EXPECT_FALSE (t8_forest_element_is_leaf_or_ghost (forest, not_leaf, itree, 0)); } } scheme->element_destroy (tree_class, 1, ¬_leaf); } } -TEST_P (element_is_leaf, element_is_leaf) +TEST_P (element_is_leaf_or_ghost, element_is_leaf) +{ + t8_test_element_is_leaf_for_forest (forest); +} + +TEST_P (element_is_leaf_or_ghost, element_is_leaf_adapt) +{ + t8_test_element_is_leaf_for_forest (forest_adapt); +} + +void +t8_test_element_is_ghost_for_forest (t8_forest_t forest) +{ + const t8_locidx_t num_ghost_trees = t8_forest_get_num_ghost_trees (forest); + const t8_scheme *scheme = t8_forest_get_scheme (forest); + for (t8_locidx_t ighost_tree = 0; ighost_tree < num_ghost_trees; ++ighost_tree) { + const t8_locidx_t num_elements_in_tree = t8_forest_ghost_tree_num_leaf_elements (forest, ighost_tree); + const t8_eclass_t tree_class = t8_forest_get_tree_class (forest, ighost_tree); + /* Allocate memory to build a non-ghost element. */ + t8_element_t *not_ghost; + scheme->element_new (tree_class, 1, ¬_ghost); + /* Iterate over all the tree's ghost elements, check whether the ghost + * is correctly identified by t8_forest_element_is_ghost and t8_forest_element_is_leaf_or_ghost, + * build its parent and its first child (if they exist), and verify + * that t8_forest_element_is_leaf and t8_forest_element_is_leaf_or_ghost returns false. */ + for (t8_locidx_t ielement = 0; ielement < num_elements_in_tree; ++ielement) { + const t8_element_t *ghost_element = t8_forest_ghost_get_leaf_element (forest, ighost_tree, ielement); + EXPECT_TRUE (t8_forest_element_is_ghost (forest, ghost_element, ighost_tree)); + EXPECT_TRUE (t8_forest_element_is_leaf_or_ghost (forest, ghost_element, ighost_tree, 1)); + /* Compute parent and first child of element and check that they are not in the tree */ + const int element_level = scheme->element_get_level (tree_class, ghost_element); + if (element_level > 0) { + scheme->element_get_parent (tree_class, ghost_element, not_ghost); + EXPECT_FALSE (t8_forest_element_is_ghost (forest, not_ghost, ighost_tree)); + EXPECT_FALSE (t8_forest_element_is_leaf_or_ghost (forest, not_ghost, ighost_tree, 1)); + } + if (element_level < scheme->get_maxlevel (tree_class)) { + scheme->element_get_child (tree_class, ghost_element, 0, not_ghost); + EXPECT_FALSE (t8_forest_element_is_ghost (forest, not_ghost, ighost_tree)); + EXPECT_FALSE (t8_forest_element_is_leaf_or_ghost (forest, not_ghost, ighost_tree, 1)); + } + } + scheme->element_destroy (tree_class, 1, ¬_ghost); + } +} + +TEST_P (element_is_leaf_or_ghost, element_is_ghost) { t8_test_element_is_leaf_for_forest (forest); } -TEST_P (element_is_leaf, element_is_leaf_adapt) +TEST_P (element_is_leaf_or_ghost, element_is_ghost_adapt) { t8_test_element_is_leaf_for_forest (forest_adapt); } -TEST_P (element_is_leaf_hybrid, element_is_leaf) +TEST_P (element_is_leaf_or_ghost_hybrid, element_is_leaf) { t8_test_element_is_leaf_for_forest (forest); } -TEST_P (element_is_leaf_hybrid, element_is_leaf_adapt) +TEST_P (element_is_leaf_or_ghost_hybrid, element_is_leaf_adapt) { t8_test_element_is_leaf_for_forest (forest_adapt); } -/* Define a lambda to beatify gtest output for tuples . +/* Define a lambda to beautify gtest output for tuples . * This will set the correct level and cmesh name as part of the test case name. */ auto pretty_print_eclass_scheme_and_level = [] (const testing::TestParamInfo, int>> &info) { @@ -204,8 +228,9 @@ auto pretty_print_eclass_scheme_and_level return scheme + "_" + eclass + level; }; -INSTANTIATE_TEST_SUITE_P (t8_gtest_element_is_leaf, element_is_leaf, +INSTANTIATE_TEST_SUITE_P (t8_gtest_element_is_leaf_or_ghost_hybrid, element_is_leaf_or_ghost_hybrid, + AllSchemeCollections, print_scheme); + +INSTANTIATE_TEST_SUITE_P (t8_gtest_element_is_leaf_or_ghost, element_is_leaf_or_ghost, testing::Combine (AllSchemes, testing::Range (0, T8_IS_LEAF_MAX_LVL)), pretty_print_eclass_scheme_and_level); - -INSTANTIATE_TEST_SUITE_P (t8_gtest_element_is_leaf_hybrid, element_is_leaf_hybrid, AllSchemeCollections, print_scheme); diff --git a/test/t8_forest/t8_gtest_forest_face_normal.cxx b/test/t8_forest/t8_gtest_forest_face_normal.cxx index 13f6d1a4a8..8b3082c027 100644 --- a/test/t8_forest/t8_gtest_forest_face_normal.cxx +++ b/test/t8_forest/t8_gtest_forest_face_normal.cxx @@ -82,16 +82,15 @@ TEST_P (class_forest_face_normal, back_and_forth) /* Get all face neighbors */ - t8_element_t **neighbors; + const t8_element_t **neighbors; int num_neighbors; - const int forest_is_balanced = 1; t8_eclass_t neigh_eclass; int *dual_faces; t8_locidx_t *neigh_ids; t8_gloidx_t gneightree; t8_forest_leaf_face_neighbors_ext (forest, itree, element, &neighbors, iface, &dual_faces, &num_neighbors, - &neigh_ids, &neigh_eclass, forest_is_balanced, &gneightree, NULL); + &neigh_ids, &neigh_eclass, &gneightree, NULL); /* Iterate and compute their facenormal */ for (int ineigh = 0; ineigh < num_neighbors; ineigh++) { @@ -104,7 +103,6 @@ TEST_P (class_forest_face_normal, back_and_forth) } 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_half_neighbors.cxx b/test/t8_forest/t8_gtest_half_neighbors.cxx index 835c0229f3..e90bad8c5a 100644 --- a/test/t8_forest/t8_gtest_half_neighbors.cxx +++ b/test/t8_forest/t8_gtest_half_neighbors.cxx @@ -96,19 +96,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); @@ -120,12 +119,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..92f626fb9a --- /dev/null +++ b/test/t8_forest/t8_gtest_leaf_face_neighbors.cxx @@ -0,0 +1,354 @@ +/* + This file is part of t8code. + t8code is a C library to manage a collection (a forest) of multiple + connected adaptive space-trees of general element classes in parallel. + + Copyright (C) 2024 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "test/t8_cmesh_generator/t8_cmesh_example_sets.hxx" +#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 level = 1; + const int adapt_levels = 2; + const int max_adapt_level = level + adapt_levels; + const bool do_ghost = true; + const bool do_recursive_adapt = true; + forests[0] = t8_forest_new_uniform (cmesh, scheme, 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 }; +}; + +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_get_tree_leaf_element_array (forest, itree) + : t8_forest_ghost_get_tree_leaf_elements (forest, ghost_tree_id); + 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. + 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 + t8_locidx_t neigh_ltreeid_from_index; + // Check that neighbor index correctly yields neighbor element. + if (neigh_index < num_local_elements) { + const t8_element_t *neighbor_from_index + = t8_forest_get_leaf_element (forest, neigh_index, &neigh_ltreeid_from_index); + EXPECT_TRUE (scheme->element_is_equal (neigh_class, neighbor_from_index, neighbor)); + } + // TODO: Check neighbor index if the element is a ghost element + + // Compute the local tree id of the neighbors tree depending on whether + // it is a local tree or a ghost 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; + if (neigh_index < num_local_elements) { + EXPECT_EQ (neigh_ltreeid, neigh_ltreeid_from_index); + } // TODO: Check neighbor ltreeid if ghost tree + // 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 + fflush (stdout); + // 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); + + if (element_index < num_local_elements) { + const t8_element_t *element_from_index = t8_forest_get_leaf_element (forest, element_index, NULL); + EXPECT_EQ (element_from_index, element) + << "Neighbor neighbor element at index " << element_index << " is not original element."; + } + // TODO: Check element index if original element is a ghost element + + // 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); diff --git a/test/t8_gtest_adapt_callbacks.cxx b/test/t8_gtest_adapt_callbacks.cxx new file mode 100644 index 0000000000..bab9f26a08 --- /dev/null +++ b/test/t8_gtest_adapt_callbacks.cxx @@ -0,0 +1,58 @@ +/* + 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) 2025 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** \file t8_gtest_adapt_callbacks.cxx +* Provide forest adapt callback functions that we use in our tests. +*/ + +#include + +/* Adapt a forest such that always the first child of a + * family is refined and no other elements. This results in a highly + * imbalanced forest. + * + * This adapt callbacks requires an integer as forest user data. + * This integer is the maximum refinement level. + */ +int +t8_test_adapt_first_child (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]))); + + int level = scheme->element_get_level (eclass, elements[0]); + + /* we set a maximum refinement level as forest user data */ + int maxlevel = *(int *) t8_forest_get_user_data (forest); + if (level >= maxlevel) { + /* Do not refine after the maxlevel */ + return 0; + } + int child_id = scheme->element_get_child_id (eclass, elements[0]); + if (child_id == 1) { + return 1; + } + return 0; +} diff --git a/test/t8_gtest_adapt_callbacks.hxx b/test/t8_gtest_adapt_callbacks.hxx new file mode 100644 index 0000000000..584015a3f1 --- /dev/null +++ b/test/t8_gtest_adapt_callbacks.hxx @@ -0,0 +1,45 @@ +/* + 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) 2025 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** \file t8_gtest_adapt_callbacks.hxx +* Provide forest adapt callback functions that we use in our tests. +*/ + +#ifndef T8_GTEST_ADAPT_CALLBACKS +#define T8_GTEST_ADAPT_CALLBACKS + +#include +#include + +/* Adapt a forest such that always the first child of a + * family is refined and no other elements. This results in a highly + * imbalanced forest. + * + * This adapt callbacks requires an integer as forest user data. + * This integer is the maximum refinement level. + */ +int +t8_test_adapt_first_child (t8_forest_t forest, t8_forest_t forest_from, t8_locidx_t which_tree, + const t8_eclass_t eclass, t8_locidx_t lelement_id, const t8_scheme *scheme, + const int is_family, const int num_elements, t8_element_t *elements[]); + +#endif /* T8_GTEST_ADAPT_CALLBACKS */ diff --git a/test/t8_schemes/t8_gtest_is_ancestor.cxx b/test/t8_schemes/t8_gtest_is_ancestor.cxx new file mode 100644 index 0000000000..3f12cd7bf8 --- /dev/null +++ b/test/t8_schemes/t8_gtest_is_ancestor.cxx @@ -0,0 +1,94 @@ +/* + 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) 2025 the developers + + t8code is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + t8code is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with t8code; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +*/ + +/** \file t8_gtest_is_ancestor.cxx + * This test checks the element_is_ancestor function of the scheme interface. + * Starting from an element we build up all its ancestor up to root level and + * test whether the element_is_ancestor function returns true. + * We also test sone cases where the function returns false, namely when putting + * in an element and an ancestor in reverse order. + * + * This test is modified from the ancestor_id test. + */ + +#include +#include +#include +#include +#include "t8_gtest_dfs_base.hxx" +#include + +class class_is_ancestor: public TestDFS { + void + check_element () override + { + t8_element_t *ancestor; + scheme->element_new (eclass, 1, &ancestor); + /* Get level of current element */ + const int level = scheme->element_get_level (eclass, element); + + /* Iterate over all levels above the current element and check + if ancestor id corresponds with the child id of elem, parent, grandparent, ... */ + for (int levels_above_elem = 0; levels_above_elem < level; levels_above_elem++) { + const int ancestor_level = level - levels_above_elem; + /* Compute the ancestor by iteratively computing the parent */ + scheme->element_copy (eclass, element, ancestor); + for (int level_diff = 0; level_diff < levels_above_elem; level_diff++) { + scheme->element_get_parent (eclass, ancestor, ancestor); + } + // Check whether element_is_ancestor correctly identifies our candidate as an ancestor + EXPECT_TRUE (scheme->element_is_ancestor (eclass, ancestor, element)); + // We check that element_is_ancestor properly returns false by + // checking that element is not an ancestor of our candidate if their levels do not agree. + if (ancestor_level != level) { + EXPECT_FALSE (scheme->element_is_ancestor (eclass, element, ancestor)); + } + } + scheme->element_destroy (eclass, 1, &ancestor); + } + + protected: + void + SetUp () override + { + /* Setup DFS test */ + dfs_test_setup (); + } + void + TearDown () override + { + /* Destroy DFS test */ + dfs_test_teardown (); + } +}; + +TEST_P (class_is_ancestor, t8_recursive_dfs_is_ancestor) +{ +#if T8CODE_TEST_LEVEL >= 1 + const int maxlvl = 4; +#else + const int maxlvl = 6; +#endif + check_recursive_dfs_to_max_lvl (maxlvl); +} + +INSTANTIATE_TEST_SUITE_P (t8_gtest_is_ancestor, class_is_ancestor, AllSchemes, print_all_schemes); diff --git a/test/t8_schemes/t8_gtest_nca.cxx b/test/t8_schemes/t8_gtest_nca.cxx index 252911581a..a8f69d117f 100644 --- a/test/t8_schemes/t8_gtest_nca.cxx +++ b/test/t8_schemes/t8_gtest_nca.cxx @@ -81,6 +81,10 @@ TEST_P (nca, nca_check_shallow) scheme->element_get_nca (tree_class, desc_a, desc_b, check); /*expect equality */ EXPECT_ELEM_EQ (scheme, tree_class, check, correct_nca); + // Check against element_is_ancestor. This adds another test layer + // to both element_get_nca and element_is_ancestor. + EXPECT_TRUE (scheme->element_is_ancestor (tree_class, check, desc_a)); + EXPECT_TRUE (scheme->element_is_ancestor (tree_class, check, desc_b)); } } } @@ -125,6 +129,10 @@ TEST_P (nca, nca_check_deep) /* Expect equality of correct_nca and check for every other class */ EXPECT_ELEM_EQ (scheme, tree_class, correct_nca, check); } + // Check against element_is_ancestor. This adds another test layer + // to both element_get_nca and element_is_ancestor. + EXPECT_TRUE (scheme->element_is_ancestor (tree_class, check, desc_a)); + EXPECT_TRUE (scheme->element_is_ancestor (tree_class, check, desc_b)); } } } @@ -179,6 +187,10 @@ t8_recursive_nca_check (t8_element_t *check_nca, t8_element_t *desc_a, t8_elemen for (j = 0; j < num_children_b; j++) { scheme->element_get_child (tree_class, parent_b, j, desc_b); scheme->element_get_nca (tree_class, desc_a, desc_b, check); + // Check against element_is_ancestor. This adds another test layer + // to both element_get_nca and element_is_ancestor. + EXPECT_TRUE (scheme->element_is_ancestor (tree_class, check, desc_a)); + EXPECT_TRUE (scheme->element_is_ancestor (tree_class, check, desc_b)); if (!scheme->element_is_equal (tree_class, check_nca, check)) { level_a = scheme->element_get_level (tree_class, desc_a); @@ -307,6 +319,10 @@ TEST_P (nca, recursive_check_higher_level) scheme->element_get_nca (tree_class, parent_a, parent_b, check); EXPECT_ELEM_EQ (scheme, tree_class, parent_a, check); EXPECT_ELEM_EQ (scheme, tree_class, parent_b, check); + // Check against element_is_ancestor. This adds another test layer + // to both element_get_nca and element_is_ancestor. + EXPECT_TRUE (scheme->element_is_ancestor (tree_class, check, parent_a)); + EXPECT_TRUE (scheme->element_is_ancestor (tree_class, check, parent_b)); } } } 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);