Feature: Ancestor search

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

src/t8_forest/t8_forest_private.h

@@ -200,7 +200,7 @@ t8_forest_get_tree_leaf_element_array_mutable (const t8_forest_t forest, t8_loci
 
 /** Search for a linear element id 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.
+ * such that the element at position i has a smaller or equal 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).
@@ -213,6 +213,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 or equal 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.

src/t8_schemes/t8_default/t8_default_common/t8_default_common.hxx

@@ -191,6 +191,40 @@ struct t8_default_scheme_common: public t8_scheme_helpers<TEclass, TUnderlyingEc
     return t8_eclass_max_num_children[TEclass];
   }
 
+  /** 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] 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.
+  */
+  inline bool
+  element_is_ancestor (const t8_element_t *element_A, const t8_element_t *element_B) const
+  {
+    /* A is ancestor of B if and only if it has smaller or equal level and
+      restricted to A's level, B has the same id as A.
+
+        level(A) <= level(B) and ID(A,level(A)) == ID(B,level(B))
+    */
+    T8_ASSERT (this->underlying ().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.

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.

src/t8_schemes/t8_standalone/t8_standalone_implementation.hxx

@@ -608,6 +608,54 @@ struct t8_standalone_scheme: public t8_scheme_helpers<TEclass, t8_standalone_sch
     return 1;
   }
 
+  // Note to devs: element_is_ancestor currently cannot be static
+  //               since it uses the non-static function element_new
+  /** 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] 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_element_t *element_A, const t8_element_t *element_B) const noexcept
+  {
+    /* We compute whether A is an ancestor of B by
+
+     - If level(A) > 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<TEclass> *el_B = (const t8_standalone_element<TEclass> *) 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<TEclass> *ancestor_casted = (t8_standalone_element<TEclass> *) 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.

test/CMakeLists.txt

@@ -135,6 +135,7 @@ add_t8_cpp_test( NAME t8_gtest_cmesh_bounding_box_serial                    SOUR
 add_t8_cpp_test( NAME t8_gtest_shmem_parallel  SOURCES t8_data/t8_gtest_shmem.cxx )
 add_t8_cpp_test( NAME t8_gtest_data_pack_parallel       SOURCES t8_data/t8_gtest_data_handler.cxx t8_data/t8_data_handler_specs.cxx)
 
+add_t8_cpp_test( NAME t8_gtest_bin_search_parallel                SOURCES t8_forest/t8_gtest_bin_search.cxx t8_gtest_adapt_callbacks.cxx )
 add_t8_cpp_test( NAME t8_gtest_element_volume_serial              SOURCES t8_forest/t8_gtest_element_volume.cxx )
 add_t8_cpp_test( NAME t8_gtest_search_parallel                    SOURCES t8_forest/t8_gtest_search.cxx )
 add_t8_cpp_test( NAME t8_gtest_half_neighbors_parallel            SOURCES t8_forest/t8_gtest_half_neighbors.cxx )
@@ -147,7 +148,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 )
@@ -182,6 +183,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 )