4 files changed, 167 insertions, 0 deletions
diff --git a/Source/Utils/Algorithms/IsIn.H b/Source/Utils/Algorithms/IsIn.H
new file mode 100644
index 000000000..c9d2f477e
--- /dev/null
+++ b/Source/Utils/Algorithms/IsIn.H
@@ -0,0 +1,58 @@
+/* Copyright 2022 Andrew Myers, Luca Fedeli, Maxence Thevenet
+ * Revathi Jambunathan
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_UTILS_ALGORITHMS_ISIN_H_
+#define WARPX_UTILS_ALGORITHMS_ISIN_H_
+
+#include <algorithm>
+#include <vector>
+
+namespace utils::algorithms
+{
+    /** \brief Returns true if an item of type TE is in a vector
+     * of TV objects (provided that TE can be converted into TV), false otherwise
+     *
+     * @tparam TV the typename of the vector elements
+     * @tparam TE the typename of the item
+     *
+     * @param vect a vector of TV objects
+     * @param elem an object of type TE
+     *
+     * @return true if elem is in vect, false otherwise
+     */
+    template <typename TV, typename TE,
+        class = typename std::enable_if<std::is_convertible<TE,TV>::value>::type>
+    bool is_in(const std::vector<TV>& vect,
+               const TE& elem)
+    {
+        return (std::find(vect.begin(), vect.end(), elem) != vect.end());
+    }
+
+
+    /** \brief Returns true if any of the items of a vector<TE> is contained
+     * in another vector<TV> (provided that TE can be converted into TV)
+     *
+     * @tparam TV the typename of the first vector elements
+     * @tparam TV the typename of the second vector elements
+     *
+     * @param vect a vector of TV objects
+     * @param elems a vector of TE objects
+     *
+     * @return true if any element of elems is in vect, false otherwise
+     */
+    template <typename TV, typename TE,
+        class = typename std::enable_if<std::is_convertible<TE,TV>::value>::type>
+    bool any_of_is_in(const std::vector<TV>& vect,
+               const std::vector<TE>& elems)
+    {
+        return std::any_of(elems.begin(), elems.end(),
+            [&](const auto elem){return is_in(vect, elem);});
+    }
+}
+
+#endif //WARPX_UTILS_ALGORITHMS_ISIN_H_
diff --git a/Source/Utils/Algorithms/LinearInterpolation.H b/Source/Utils/Algorithms/LinearInterpolation.H
new file mode 100644
index 000000000..32fdf7a6e
--- /dev/null
+++ b/Source/Utils/Algorithms/LinearInterpolation.H
@@ -0,0 +1,59 @@
+/* Copyright 2022 Luca Fedeli
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_UTILS_ALGORITHMS_LINEAR_INTERPOLATION_H_
+#define WARPX_UTILS_ALGORITHMS_LINEAR_INTERPOLATION_H_
+
+#include <AMReX_Extension.H>
+#include <AMReX_GpuQualifiers.H>
+
+namespace utils::algorithms
+{
+    /** \brief Performs a linear interpolation
+     *
+     * Performs a linear interpolation at x given the 2 points
+     * (x0, f0) and (x1, f1)
+     */
+    template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+    T linear_interp(T x0, T x1, T f0, T f1, T x)
+    {
+        return ((x1-x)*f0 + (x-x0)*f1)/(x1-x0);
+    }
+
+    /** \brief Performs a bilinear interpolation
+     *
+     * Performs a bilinear interpolation at (x,y) given the 4 points
+     * (x0, y0, f00), (x0, y1, f01), (x1, y0, f10), (x1, y1, f11).
+     */
+    template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+    T bilinear_interp(T x0, T x1, T y0, T y1, T f00, T f01, T f10, T f11, T x, T y)
+    {
+        const T fx0 = linear_interp(x0, x1, f00, f10, x);
+        const T fx1 = linear_interp(x0, x1, f01, f11, x);
+        return linear_interp(y0, y1, fx0, fx1, y);
+    }
+
+    /** \brief Performs a trilinear interpolation
+     *
+     * Performs a trilinear interpolation at (x,y,z) given the 8 points
+     * (x0, y0, z0, f000), (x0, y0, z1, f001), (x0, y1, z0, f010), (x0, y1, z1, f011),
+     * (x1, y0, z0, f100), (x1, y0, z1, f101), (x1, y1, z0, f110), (x1, y1, z1, f111)
+     */
+    template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+    T trilinear_interp(T x0, T x1,T y0, T y1, T z0, T z1,
+        T f000, T f001, T f010, T f011, T f100, T f101, T f110, T f111,
+        T x, T y, T z)
+    {
+        const T fxy0 = bilinear_interp(
+            x0, x1, y0, y1, f000, f010, f100, f110, x, y);
+        const T fxy1 = bilinear_interp(
+            x0, x1, y0, y1, f001, f011, f101, f111, x, y);
+        return linear_interp(z0, z1, fxy0, fxy1, z);
+    }
+}
+
+#endif //WARPX_UTILS_ALGORITHMS_LINEAR_INTERPOLATION_H_
diff --git a/Source/Utils/Algorithms/Make.package b/Source/Utils/Algorithms/Make.package
new file mode 100644
index 000000000..04638494d
--- /dev/null
+++ b/Source/Utils/Algorithms/Make.package
@@ -0,0 +1 @@
+VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Utils/Algorithms
diff --git a/Source/Utils/Algorithms/UpperBound.H b/Source/Utils/Algorithms/UpperBound.H
new file mode 100644
index 000000000..8a528971a
--- /dev/null
+++ b/Source/Utils/Algorithms/UpperBound.H
@@ -0,0 +1,49 @@
+/* Copyright 2022 Luca Fedeli
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_UTILS_ALGORITHMS_UPPER_BOUND_H_
+#define WARPX_UTILS_ALGORITHMS_UPPER_BOUND_H_
+
+#include <AMReX_Extension.H>
+#include <AMReX_GpuQualifiers.H>
+
+namespace utils::algorithms
+{
+
+    /** \brief Returns a pointer to the first element in the range [first, last) that is greater than val
+     *
+     * A re-implementation of the upper_bound algorithm suitable for GPU kernels.
+     *
+     * @param first: pointer to left limit of the range to consider
+     * @param last: pointer to right limit of the range to consider
+     * @param val: value to compare the elements of [first, last) to
+     */
+    template<typename T> AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+    const T* upper_bound(const T* first, const T* last, const T& val)
+    {
+        const T* it;
+        size_t count, step;
+        count = last-first;
+        while(count>0){
+            it = first;
+            step = count/2;
+            it += step;
+            if (!(val<*it)){
+                first = ++it;
+                count -= step + 1;
+            }
+            else{
+                count = step;
+            }
+        }
+
+        return first;
+    }
+
+}
+
+#endif //WARPX_UTILS_ALGORITHMS_UPPER_BOUND_H_