Feature: extend leaf face neighbor to non-balanced forests and ghosts [N/N]

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);
               }
 

example/forest/t8_test_face_iterate.cxx

@@ -28,6 +28,7 @@
 #include <t8_forest/t8_forest_general.h>
 #include <t8_forest/t8_forest_io.h>
 #include <t8_forest/t8_forest_geometrical.h>
+#include <t8_forest/t8_forest_ghost.h>
 #include <t8_cmesh/t8_cmesh.h>
 #include <t8_cmesh/t8_cmesh_io/t8_cmesh_readmshfile.h>
 #include <t8_cmesh/t8_cmesh_examples.h>
@@ -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) {

src/t8_data/t8_containers.cxx

@@ -27,6 +27,7 @@
 #include <sc_containers.h>
 #include <t8_data/t8_containers.h>
 #include <t8_schemes/t8_scheme.hxx>
+#include <t8_forest/t8_forest_private.h>
 
 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)
 {

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,