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

NEWS.md

@@ -21,6 +21,24 @@ No! If your code is only a couple of lines long AND has very little impact on th
 
 # User Updates for the upcoming t8code release (February 2026 - version format unclear)
 
+## Updates to t8_forest_leaf_face_neighbors
+
+The t8_forest_leaf_face_neighbors and t8_forest_leaf_face_neighbors_ext functions (see t8_forest_general.h) have undergone a major update.
+
+- The forest does not need to be balanced anymore to compute the leaf face neighbors. Thus, there can be any arbitrary number of leaf face neighbors for a given element.
+The last argument (forest_is_balanced) was removed.
+
+- The allocation behaviour changed. element_destroy should not be called anymore.
+The function now returns pointers to t8code internal elements in the neighbor_leaves array.
+To free the memory use
+```C++
+if (num_neighbors > 0) {
+   T8_FREE (pneighbor_leaves);
+   T8_FREE (pelement_indices);
+   T8_FREE (dual_faces);
+}
+```
+
 ## Using structs instead of classes
 
 We decided to use `struct` instead of `class` throughout the entire codebase.

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_ghost_get_tree_leaf_elements (forest, ghost_tree_id)
+                                                : t8_forest_tree_get_leaf_elements (forest, itree);
+    // Get the first and last element
+    const t8_element_t *first_el = (const t8_element_t *) t8_element_array_index_locidx (leaf_elements, 0);
     const t8_element_t *last_el
-      = t8_forest_get_leaf_element_in_tree (forest, itree, t8_forest_get_tree_num_leaf_elements (forest, itree) - 1);
+      = (const t8_element_t *) t8_element_array_index_locidx (leaf_elements, num_tree_elements - 1);
+
+    t8_element_t *nca;
     scheme->element_new (eclass, 1, &nca);
     scheme->element_get_nca (eclass, first_el, last_el, nca);
-    leaf_elements = t8_forest_tree_get_leaf_elements (forest, itree);
 
-    for (iface = 0; iface < scheme->element_get_num_faces (eclass, nca); iface++) {
+    //
+    for (int iface = 0; iface < scheme->element_get_num_faces (eclass, nca); iface++) {
       udata.count = 0;
-      t8_forest_iterate_faces (forest, itree, nca, iface, leaf_elements, &udata, 0, t8_test_fiterate_callback);
+      t8_forest_iterate_faces (forest, itree, nca, iface, leaf_elements, 0, t8_test_fiterate_callback, &udata);
       t8_debugf ("Leaf elements at face %i:\t%i\n", iface, udata.count);
     }
     scheme->element_destroy (eclass, 1, &nca);
@@ -139,6 +150,7 @@ t8_test_fiterate_refine_and_partition (t8_cmesh_t cmesh, int level, sc_MPI_Comm
   /* partition the adapted forest */
   t8_forest_init (&forest_partition);
   t8_forest_set_partition (forest_partition, forest_adapt, 0);
+  t8_forest_set_ghost (forest_partition, 1, T8_GHOST_FACES);
   t8_forest_commit (forest_partition);
   t8_debugf ("Created ghost structure with %li ghost elements.\n", (long) t8_forest_get_num_ghosts (forest_partition));
   if (!no_vtk) {

mesh_handle/element.hxx

@@ -345,15 +345,14 @@ class element: public TCompetences<element<TMeshClass, TCompetences...>>... {
         return this->m_neighbors[face];
       }
     }
-
-    t8_element_t** neighbors; /**< Neighboring elements. */
-    int* dual_faces_internal; /**< Face indices of the neighbor elements. */
-    int num_neighbors;        /**< Number of neighboring elements. */
-    t8_locidx_t* neighids;    /**< Neighboring elements ids. */
-    t8_eclass_t neigh_class;  /**< Neighboring elements tree class. */
+    const t8_element_t** neighbors; /**< Neighboring elements. */
+    int* dual_faces_internal;       /**< Face indices of the neighbor elements. */
+    int num_neighbors;              /**< Number of neighboring elements. */
+    t8_locidx_t* neighids;          /**< Neighboring elements ids. */
+    t8_eclass_t neigh_class;        /**< Neighboring elements tree class. */
 
     t8_forest_leaf_face_neighbors (m_mesh->m_forest, m_tree_id, m_element, &neighbors, face, &dual_faces_internal,
-                                   &num_neighbors, &neighids, &neigh_class, t8_forest_is_balanced (m_mesh->m_forest));
+                                   &num_neighbors, &neighids, &neigh_class);
     if (dual_faces) {
       dual_faces->get ().assign (dual_faces_internal, dual_faces_internal + num_neighbors);
     }
@@ -369,7 +368,6 @@ class element: public TCompetences<element<TMeshClass, 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);