Merge pull request #1162 from DLR-AMR/enhancement-use_std_lower_bound

Enhancement #1089: Replace the implemented binary search by a binary search from the standard library

src/Makefile.am

@@ -44,7 +44,8 @@ libt8_installed_headers_cmesh = \
   src/t8_cmesh/t8_cmesh_types.h \
   src/t8_cmesh/t8_cmesh_stash.h
 libt8_installed_headers_data = \
-  src/t8_data/t8_shmem.h src/t8_data/t8_containers.h
+  src/t8_data/t8_shmem.h src/t8_data/t8_containers.h \
+  src/t8_data/t8_element_array_iterator.hxx
 libt8_installed_headers_forest = \
   src/t8_forest/t8_forest.h \
   src/t8_forest/t8_forest_general.h \

src/t8_data/t8_element_array_iterator.hxx

@@ -0,0 +1,241 @@
+/*
+  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.
+*/
+
+/** \file t8_element_array_iterator.hxx
+ * We define a Random-Access-Iterator wrapping around a \a t8_element_array_t that stores pointers to elements of a given
+ * eclass scheme.
+ */
+
+#ifndef T8_ELEMENT_ARRAY_ITERATOR_HXX
+#define T8_ELEMENT_ARRAY_ITERATOR_HXX
+
+#include <t8.h>
+#include <t8_element.h>
+#include <t8_data/t8_containers.h>
+
+#include <cstddef>
+#include <iterator>
+
+/**
+ * \brief This class resembles a Random-Access-Iterator wrapper for the content of a \a t8_element_array_t.
+ * The iterator can be dereferenced in order to receive a \a t8_element_t pointer to the corresponding index.
+ * These iterators may be used in algorithms of the standard library (e.g. std::lower_bound,...) or to iterate
+ * through a \a t8_element_array_t.
+ * Since the actual data is serialized to char-bytes in the underlying array, we reconstruct a \a t8_element_t pointer
+ * and let the dereference operator return it as a value_type. Therefore, read-only operations on the
+ * \a t8_element_array_t are possible.
+ */
+class t8_element_array_iterator {
+
+ private:
+  const t8_eclass_scheme_c* scheme; /*!< The scheme of the elements residing within the array. */
+  const sc_array_t* elements;       /*!< A pointer to the actual serialized array of element pointers. */
+  t8_locidx_t current_index { 0 };  /*!< The index the iterator currently points to. */
+
+ public:
+  using iterator_category = std::random_access_iterator_tag;
+  using difference_type = std::ptrdiff_t;
+  using pointer = t8_element_t**;
+  using value_type = t8_element_t*;
+  using reference = t8_element_t*&;
+
+  /* Constructors */
+  t8_element_array_iterator () = delete;
+  t8_element_array_iterator (const t8_element_array_t* element_array, const t8_locidx_t position)
+    : scheme { t8_element_array_get_scheme (element_array) }, elements { t8_element_array_get_array (element_array) },
+      current_index { position } {};
+
+  /* Copy/Move Constructors/Assignment-Operators */
+  t8_element_array_iterator (const t8_element_array_iterator& other) = default;
+  t8_element_array_iterator&
+  operator= (const t8_element_array_iterator& other)
+    = default;
+  t8_element_array_iterator (t8_element_array_iterator&& other) = default;
+  t8_element_array_iterator&
+  operator= (t8_element_array_iterator&& other)
+    = default;
+
+  /* Destructor */
+  ~t8_element_array_iterator () = default;
+
+  /* Dereferencing operator of the iterator wrapper returning a value_type (a t8_element_t-pointer
+   * casted from the serialized char-bytes in the underlying sc_array_t). */
+  value_type
+  operator* ()
+  {
+    T8_ASSERT (current_index >= 0 && static_cast<size_t> (current_index) < elements->elem_count);
+    return static_cast<t8_element_t*> (t8_sc_array_index_locidx (elements, current_index));
+  };
+
+  value_type
+  operator[] (const difference_type n) const
+  {
+    T8_ASSERT (n >= 0 && static_cast<size_t> (n) < elements->elem_count);
+    return static_cast<t8_element_t*> (t8_sc_array_index_locidx (elements, n));
+  };
+
+  /* Pre- and Postfix increment */
+  t8_element_array_iterator&
+  operator++ ()
+  {
+    ++current_index;
+    return *this;
+  };
+  t8_element_array_iterator
+  operator++ (int)
+  {
+    t8_element_array_iterator tmp_iterator (*this);
+    ++(*this);
+    return tmp_iterator;
+  };
+
+  /* Pre- and Postfix decrement */
+  t8_element_array_iterator&
+  operator-- ()
+  {
+    --current_index;
+    return *this;
+  };
+  t8_element_array_iterator
+  operator-- (int)
+  {
+    t8_element_array_iterator tmp_iterator (*this);
+    --(*this);
+    return tmp_iterator;
+  };
+
+  /* Arithmetic assignment operators */
+  t8_element_array_iterator&
+  operator+= (const difference_type n)
+  {
+    current_index += n;
+    return *this;
+  }
+  t8_element_array_iterator&
+  operator-= (const difference_type n)
+  {
+    current_index -= n;
+    T8_ASSERT (current_index < 0);
+    return *this;
+  }
+  /* Comparison operators */
+  friend bool
+  operator== (const t8_element_array_iterator& iter1, const t8_element_array_iterator& iter2)
+  {
+    return (iter1.elements->array == iter2.elements->array && iter1.current_index == iter2.current_index);
+  };
+  friend bool
+  operator!= (const t8_element_array_iterator& iter1, const t8_element_array_iterator& iter2)
+  {
+    return (iter1.elements->array != iter2.elements->array || iter1.current_index != iter2.current_index);
+  };
+  friend bool
+  operator< (const t8_element_array_iterator& lhs, const t8_element_array_iterator& rhs)
+  {
+    T8_ASSERT (lhs.elements->array == rhs.elements->array);
+    return lhs.current_index < rhs.current_index;
+  }
+  friend bool
+  operator> (const t8_element_array_iterator& lhs, const t8_element_array_iterator& rhs)
+  {
+    T8_ASSERT (rhs.elements->array == lhs.elements->array);
+    return rhs.current_index < lhs.current_index;
+  }
+  friend bool
+  operator<= (const t8_element_array_iterator& lhs, const t8_element_array_iterator& rhs)
+  {
+    return !(rhs < lhs);
+  }
+  friend bool
+  operator>= (const t8_element_array_iterator& lhs, const t8_element_array_iterator& rhs)
+  {
+    return !(lhs < rhs);
+  }
+
+  /* Arithmetic operators */
+  friend t8_element_array_iterator
+  operator+ (const t8_element_array_iterator& iter, const difference_type n)
+  {
+    t8_element_array_iterator tmp_iterator (iter);
+    tmp_iterator += n;
+    return tmp_iterator;
+  }
+  friend t8_element_array_iterator
+  operator+ (const difference_type n, const t8_element_array_iterator& iter)
+  {
+    return iter + n;
+  }
+  friend t8_element_array_iterator
+  operator- (const t8_element_array_iterator& iter, const difference_type n)
+  {
+    t8_element_array_iterator tmp_iterator (iter);
+    tmp_iterator -= n;
+    return tmp_iterator;
+  }
+  friend difference_type
+  operator- (const t8_element_array_iterator& lhs, const t8_element_array_iterator& rhs)
+  {
+    T8_ASSERT (lhs.elements->array == rhs.elements->array);
+    return lhs.current_index - rhs.current_index;
+  }
+
+  /* Return the index within the array the iterator currently points to [0,...,size]. */
+  t8_locidx_t
+  GetCurrentIndex () const
+  {
+    return current_index;
+  };
+
+  /* Compute the linear id at a given level for the element the iterator points to. */
+  t8_linearidx_t
+  GetLinearIDAtLevel (const int level)
+  {
+    T8_ASSERT (current_index >= 0 && static_cast<size_t> (current_index) < elements->elem_count);
+    return scheme->t8_element_get_linear_id (*(*this), level);
+  };
+};
+
+/**
+ * \brief Get an Random-Access-Iterator to the first element pointer within an element array.
+ * 
+ * \param[in] element_array The element array from which the begin()-Random-Access-Iterator should be taken.
+ * \return t8_element_array_iterator The iterator to the start of the element array.
+ */
+inline t8_element_array_iterator
+t8_element_array_begin (const t8_element_array_t* element_array)
+{
+  return t8_element_array_iterator (element_array, 0);
+}
+
+/**
+ * \brief Get an Random-Access-Iterator to the one after the last element pointer within an element array.
+ * 
+ * \param[in] element_array The element array from which the end()-Random-Access-Iterator should be taken.
+ * \return t8_element_array_iterator The iterator to the end of the element array.
+ */
+inline t8_element_array_iterator
+t8_element_array_end (const t8_element_array_t* element_array)
+{
+  return t8_element_array_iterator (element_array, t8_element_array_get_count (element_array));
+}
+
+#endif /* !T8_ELEMENT_ARRAY_ITERATOR_HXX */

src/t8_forest/t8_forest.cxx

@@ -43,6 +43,9 @@
 #include <t8_geometry/t8_geometry_implementations/t8_geometry_linear.h>
 #include <t8_geometry/t8_geometry_implementations/t8_geometry_linear_axis_aligned.h>
 #endif
+#include <t8_data/t8_element_array_iterator.hxx>
+
+#include <algorithm>
 
 /* We want to export the whole implementation to be callable from "C" */
 T8_EXTERN_C_BEGIN ();
@@ -1438,50 +1441,37 @@ t8_forest_copy_trees (t8_forest_t forest, t8_forest_t from, int copy_elements)
   }
 }
 
-/* Search for a linear element id (at forest->maxlevel) in a sorted array of
+/** \brief Search for a linear element id (at forest->maxlevel) 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.
  * If no such i exists, return -1.
  */
-/* TODO: should return t8_locidx_t */
 static t8_locidx_t
 t8_forest_bin_search_lower (const t8_element_array_t *elements, const t8_linearidx_t element_id, const int maxlevel)
 {
-  t8_linearidx_t query_id;
-  t8_locidx_t low, high, guess;
-
   const t8_eclass_scheme_c *ts = t8_element_array_get_scheme (elements);
   /* At first, we check whether any element has smaller id than the
    * given one. */
   const t8_element_t *query = t8_element_array_index_int (elements, 0);
-  query_id = ts->t8_element_get_linear_id (query, maxlevel);
+  const t8_linearidx_t query_id = ts->t8_element_get_linear_id (query, maxlevel);
   if (query_id > element_id) {
-    /* No element has id smaller than the given one */
+    /* No element has id smaller than the given one. */
     return -1;
   }
 
-  /* We now perform the binary search */
-  low = 0;
-  high = t8_element_array_get_count (elements) - 1;
-  while (low < high) {
-    guess = (low + high + 1) / 2;
-    query = t8_element_array_index_int (elements, guess);
-    query_id = ts->t8_element_get_linear_id (query, maxlevel);
-    if (query_id == element_id) {
-      /* we are done */
-      return guess;
-    }
-    else if (query_id > element_id) {
-      /* look further left */
-      high = guess - 1;
-    }
-    else {
-      /* look further right, but keep guess in the search range */
-      low = guess;
-    }
-  }
-  T8_ASSERT (low == high);
-  return low;
+  /* We search for the first element in the array that is greater than the given element id. */
+  auto elem_iter = std::upper_bound (
+    t8_element_array_begin (elements), t8_element_array_end (elements), element_id,
+    [&maxlevel, &ts] (const t8_linearidx_t element_id_, const t8_element_array_iterator::value_type &elem_ptr) {
+      return (element_id_ < ts->t8_element_get_linear_id (elem_ptr, maxlevel));
+    });
+
+  /* After we found the element with an id greater than the given one, we are able to jump one index back.
+   * This guarantees us that the element at (index - 1) is smaller or equal to the given element id.
+   * 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. In that case, we want to return the last index from the element
+   * array. */
+  return elem_iter.GetCurrentIndex () - 1;
 }
 
 t8_eclass_t