add interface for replicated partition search

doc/author_brandt.txt

@@ -0,0 +1 @@
+I place my contributions to t8code under the FreeBSD license. Hannes Brandt ([email protected])

src/CMakeLists.txt

@@ -112,7 +112,8 @@ target_sources( T8 PRIVATE
     t8_forest/t8_forest_ghost.cxx 
     t8_forest/t8_forest_iterate.cxx 
     t8_forest/t8_forest_balance.cxx 
-    t8_forest/t8_forest_netcdf.cxx 
+    t8_forest/t8_forest_netcdf.cxx
+    t8_forest/t8_forest_search/t8_forest_search.cxx 
     t8_geometry/t8_geometry.cxx 
     t8_geometry/t8_geometry_helpers.c 
     t8_geometry/t8_geometry_base.cxx 

src/Makefile.am

@@ -53,7 +53,9 @@ libt8_installed_headers_forest = \
   src/t8_forest/t8_forest_profiling.h \
   src/t8_forest/t8_forest_io.h \
   src/t8_forest/t8_forest_adapt.h \
-  src/t8_forest/t8_forest_iterate.h src/t8_forest/t8_forest_partition.h
+  src/t8_forest/t8_forest_iterate.h \
+  src/t8_forest/t8_forest_partition.h \
+  src/t8_forest/t8_forest_search/t8_forest_search.hxx
 libt8_installed_headers_geometry = \
   src/t8_geometry/t8_geometry.h \
   src/t8_geometry/t8_geometry_handler.hxx \
@@ -139,6 +141,7 @@ libt8_compiled_sources = \
   src/t8_forest/t8_forest_partition.cxx src/t8_forest/t8_forest.cxx \
   src/t8_forest/t8_forest_private.c \
   src/t8_forest/t8_forest_ghost.cxx src/t8_forest/t8_forest_iterate.cxx \
+  src/t8_forest/t8_forest_search/t8_forest_search.cxx \
   src/t8_version.c \
   src/t8_vtk.c src/t8_forest/t8_forest_balance.cxx \
   src/t8_forest/t8_forest_netcdf.cxx \

src/t8_forest/t8_forest.cxx

@@ -3719,7 +3719,6 @@ t8_eclass_t
 t8_forest_get_tree_class (const t8_forest_t forest, const t8_locidx_t ltreeid)
 {
   t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (forest);
-
   T8_ASSERT (0 <= ltreeid && ltreeid < num_local_trees + t8_forest_get_num_ghost_trees (forest));
   if (ltreeid < num_local_trees) {
     /* The id belongs to a local tree */

src/t8_forest/t8_forest_iterate.cxx

@@ -530,4 +530,11 @@ t8_forest_iterate_replace (t8_forest_t forest_new, t8_forest_t forest_old, t8_fo
   t8_global_productionf ("Done t8_forest_iterate_replace\n");
 }
 
+void
+t8_forest_search_partition (const t8_forest_t forest, t8_forest_partition_search_fn search_fn,
+                            t8_forest_partition_query_fn query_fn, sc_array_t *queries)
+{
+  T8_ASSERT (0); /* not implemented yet */
+}
+
 T8_EXTERN_C_END ();

src/t8_forest/t8_forest_iterate.h

@@ -79,6 +79,48 @@ typedef void (*t8_forest_query_fn) (t8_forest_t forest, const t8_locidx_t ltreei
                                     const int is_leaf, const t8_element_array_t *leaf_elements,
                                     const t8_locidx_t tree_leaf_index, sc_array_t *queries, sc_array_t *query_indices,
                                     int *query_matches, const size_t num_active_queries);
+/**
+ * A call-back function used by \ref t8_forest_search_partition describing a search-criterion. Is called on an element
+ * and the search criterion should be checked on that element. Return true if the search criterion is met, false
+ * otherwise.
+ *
+ * \param[in] forest              the forest
+ * \param[in] ltreeid             the local tree id of the current tree in the cmesh. Since the cmesh has to be
+ *                                replicated, it coincides with the global tree id.
+ * \param[in] element             the element for which the search criterion is checked
+ * \param[in] pfirst              the first processor that owns part of \a element. Guaranteed to be non-empty.
+ * \param[in] plast               the last processor that owns part of \a element. Guaranteed to be non-empty.
+ * \returns                       non-zero if the search criterion is met, zero otherwise.
+ */
+typedef int (*t8_forest_partition_search_fn) (const t8_forest_t forest, const t8_locidx_t ltreeid,
+                                              const t8_element_t *element, const int pfirst, const int plast);
+
+/**
+ * A call-back function used by \ref t8_forest_search_partition for queries. Is called on an element and all queries are
+ * checked on that element. All positive queries are passed further down to the children of the element. The results of
+ * the check are stored in \a query_matches.
+ *
+ * \param[in] forest              the forest
+ * \param[in] ltreeid             the local tree id of the current tree in the cmesh. Since the cmesh has to be
+ *                                replicated, it coincides with the global tree id.
+ * \param[in] element             the element for which the query is executed
+ * \param[in] pfirst              the first processor that owns part of \a element. Guaranteed to be non-empty.
+ * \param[in] plast               the last processor that owns part of \a element. Guaranteed to be non-empty.
+ *                                if this is equal to \a pfirst, then the recursion will stop for
+ *                                \a element's branch after this function returns.
+ * \param[in] queries             an array of queries that are checked by the function
+ * \param[in] query_indices       an array of size_t entries, where each entry is an index of a query in \a queries.
+ * \param[in, out] query_matches  an array of length \a num_active_queries.
+ *                                If the element is not a leaf must be set to true or false at the i-th index for
+ *                                each query, specifying whether the element 'matches' the query of the i-th query
+ *                                index or not. When the element is a leaf we can return before all entries are set.
+ * \param[in] num_active_queries  The number of currently active queries (equals the number of entries of
+ *                                \a query_matches and entries of \a query_indices).
+ */
+typedef void (*t8_forest_partition_query_fn) (const t8_forest_t forest, const t8_locidx_t ltreeid,
+                                              const t8_element_t *element, const int pfirst, const int plast,
+                                              void *queries, sc_array_t *query_indices, int *query_matches,
+                                              const size_t num_active_queries);
 
 T8_EXTERN_C_BEGIN ();
 
@@ -134,6 +176,25 @@ t8_forest_search (t8_forest_t forest, t8_forest_search_fn search_fn, t8_forest_q
 void
 t8_forest_iterate_replace (t8_forest_t forest_new, t8_forest_t forest_old, t8_forest_replace_t replace_fn);
 
+/**
+ * Perform a top-down search of the global partition, executing a callback on
+ * each intermediate element. The search will enter each tree at least once.
+ * The recursion will only go down branches that are split between multiple processors.
+ * This is not a collective function. It does not communicate.
+ * The function expects the coarse mesh to be replicated.
+ * If the callback returns false for an element, its descendants
+ * are not further searched.
+ * To pass user data to \b search_fn function use \ref t8_forest_set_user_data
+ *
+ * \param[in] forest              the forest to be searched
+ * \param[in] search_fn           a search callback function called on elements
+ * \param[in] query_fn            a query callback function called for all active queries of an element
+ * \param[in,out] queries         an array of queries that are checked by the function
+ */
+void
+t8_forest_search_partition (const t8_forest_t forest, t8_forest_partition_search_fn search_fn,
+                            t8_forest_partition_query_fn query_fn, sc_array_t *queries);
+
 T8_EXTERN_C_END ();
 
 #endif /* !T8_FOREST_ITERATE_H */

src/t8_forest/t8_forest_search/t8_forest_search.cxx

@@ -0,0 +1,142 @@
+/*
+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 "t8_forest/t8_forest_search/t8_forest_search.hxx"
+#include <t8_forest/t8_forest_iterate.h>
+#include <t8_forest/t8_forest_types.h>
+#include <t8_forest/t8_forest_general.h>
+#include <t8_element.hxx>
+
+void
+t8_search_base::search_recursion (const t8_locidx_t ltreeid, t8_element_t *element, const t8_eclass_scheme_c *ts,
+                                  t8_element_array_t *leaf_elements, const t8_locidx_t tree_lindex_of_first_leaf)
+{
+  /* Assertions to check for necessary requirements */
+  /* The forest must be committed */
+  T8_ASSERT (t8_forest_is_committed (this->forest));
+  /* The tree must be local */
+  T8_ASSERT (0 <= ltreeid && ltreeid < t8_forest_get_num_local_trees (this->forest));
+
+  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;
+  }
+
+  if (this->stop_due_to_queries ()) {
+    return;
+  }
+
+  bool is_leaf = false;
+  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 */
+    const t8_element_t *leaf = t8_element_array_index_locidx (leaf_elements, 0);
+
+    SC_CHECK_ABORT (ts->t8_element_level (element) <= ts->t8_element_level (leaf),
+                    "Search: element level greater than leaf level\n");
+    if (ts->t8_element_level (element) == ts->t8_element_level (leaf)) {
+      T8_ASSERT (t8_forest_element_is_leaf (this->forest, leaf, ltreeid));
+      T8_ASSERT (ts->t8_element_equal (element, leaf));
+      /* The element is the leaf */
+      is_leaf = true;
+    }
+  }
+  /* Call the callback function for the element */
+  const bool ret = check_element (ltreeid, element, is_leaf, leaf_elements, tree_lindex_of_first_leaf);
+
+  if (!ret) {
+    /* The function returned false. We abort the recursion */
+    return;
+  }
+  std::vector<size_t> new_active_queries;
+  this->check_queries (new_active_queries, ltreeid, element, is_leaf, leaf_elements, tree_lindex_of_first_leaf);
+
+  if (is_leaf) {
+    return;
+  }
+
+  /* Enter the recursion (the element is definitely not a leaf at this point) */
+  /* We compute all children of E, compute their leaf arrays and call search_recursion */
+  /* allocate the memory to store the children */
+  const int num_children = ts->t8_element_num_children (element);
+  t8_element_t **children = T8_ALLOC (t8_element_t *, num_children);
+  ts->t8_element_new (num_children, children);
+  /* Memory for the indices that split the leaf_elements array */
+  size_t *split_offsets = T8_ALLOC (size_t, num_children + 1);
+  /* Compute the children */
+  ts->t8_element_children (element, num_children, children);
+  /* Split the leaves array in portions belonging to the children of element */
+  t8_forest_split_array (element, leaf_elements, split_offsets);
+  for (int ichild = 0; ichild < num_children; ichild++) {
+    /* Check if there are any leaf elements for this child */
+    const size_t indexa = split_offsets[ichild];     /* first leaf of this child */
+    const size_t indexb = split_offsets[ichild + 1]; /* first leaf of next child */
+    if (indexa < indexb) {
+      t8_element_array_t child_leaves;
+      /* There exist leaves of this child in leaf_elements,
+       * we construct an array of these leaves */
+      t8_element_array_init_view (&child_leaves, leaf_elements, indexa, indexb - indexa);
+      /* Enter the recursion */
+      search_recursion (ltreeid, children[ichild], ts, &child_leaves, indexa + tree_lindex_of_first_leaf);
+      update_queries (new_active_queries);
+    }
+  }
+
+  /* clean-up */
+  ts->t8_element_destroy (num_children, children);
+  T8_FREE (children);
+  T8_FREE (split_offsets);
+}
+
+void
+t8_search_base::search_tree (const t8_locidx_t ltreeid)
+{
+  const t8_eclass_t eclass = t8_forest_get_eclass (this->forest, ltreeid);
+  const t8_eclass_scheme_c *ts = t8_forest_get_eclass_scheme (this->forest, eclass);
+  t8_element_array_t *leaf_elements = t8_forest_tree_get_leaves (this->forest, ltreeid);
+
+  /* assert for empty tree */
+  T8_ASSERT (t8_element_array_get_count (leaf_elements) >= 0);
+  /* Get the first and last leaf of this tree */
+  const t8_element_t *first_el = t8_element_array_index_locidx (leaf_elements, 0);
+  const t8_element_t *last_el
+    = t8_element_array_index_locidx (leaf_elements, t8_element_array_get_count (leaf_elements) - 1);
+  /* Compute their nearest common ancestor */
+  t8_element_t *nca;
+  ts->t8_element_new (1, &nca);
+  ts->t8_element_nca (first_el, last_el, nca);
+
+  /* Start the top-down search */
+  this->search_recursion (ltreeid, nca, ts, leaf_elements, 0);
+
+  ts->t8_element_destroy (1, &nca);
+}
+
+void
+t8_search_base::do_search ()
+{
+  T8_ASSERT (t8_forest_is_committed (forest));
+  const t8_locidx_t num_local_trees = t8_forest_get_num_local_trees (this->forest);
+  for (t8_locidx_t itree = 0; itree < num_local_trees; itree++) {
+    this->search_tree (itree);
+  }
+}