From 58fa042265d143fa968683f0a8e3195651ecad8d Mon Sep 17 00:00:00 2001
From: Mayukh Kundu <mzk0148@auburn.edu>
Date: Fri, 4 Apr 2025 14:11:29 -0500
Subject: [PATCH 1/6] Extend center, lower_bound, upper_bound, shape, step and
 bin to higher dimensions

---
 include/flyft/mesh.h | 40 ++++++++++++------------
 src/mesh.cc          | 73 +++++++++++++++++++++++++++++++-------------
 2 files changed, 72 insertions(+), 41 deletions(-)

diff --git a/include/flyft/mesh.h b/include/flyft/mesh.h
index ab85cc0..ba70aff 100644
--- a/include/flyft/mesh.h
+++ b/include/flyft/mesh.h
@@ -19,22 +19,22 @@ class Mesh
          BoundaryType lower_bc,
          BoundaryType upper_bc);
 
-    std::shared_ptr<Mesh> slice(int start, int end) const;
+    std::shared_ptr<Mesh> slice(std::vector<int> start, std::vector<int> end) const;
 
     //! Get position on the mesh, defined as center of bin
-    double center(int i) const;
+    std::vector<double> center(std::vector<int> i) const;
 
     //! Lower bound of entire mesh
-    double lower_bound() const;
+    std::vector<double> lower_bound() const;
 
     //! Get lower bound of bin
-    double lower_bound(int i) const;
+    std::vector<double> lower_bound(std::vector<int> i) const;
 
     //! Upper bound of entire mesh
-    double upper_bound() const;
+    std::vector<double> upper_bound() const;
 
     //! Get upper bound of bin
-    double upper_bound(int i) const;
+    std::vector<double> upper_bound(int i) const;
 
     //! Get surface area of lower edge of bin
     virtual double area(int i) const = 0;
@@ -46,16 +46,16 @@ class Mesh
     virtual double volume(int i) const = 0;
 
     //! Get the bin for a coordinate
-    int bin(double x) const;
+    std::vector<int> bin(double x) const;
 
     //! Length of the mesh
     double L() const;
 
     //! Shape of the mesh
-    int shape() const;
+    std::vector<int> shape() const;
 
     //! Step size of the mesh
-    double step() const;
+    std::vector<double> step() const;
 
     //! Boundary condition on lower bound of mesh
     BoundaryType lower_boundary_condition() const;
@@ -67,31 +67,31 @@ class Mesh
 
     double asLength(int shape) const;
 
-    double integrateSurface(int idx, double j_lo, double j_hi) const;
-    double integrateSurface(int idx, const DataView<double>& j) const;
-    double integrateSurface(int idx, const DataView<const double>& j) const;
+    double integrateSurface(std::vector<int> idx, double j_lo, double j_hi) const;
+    double integrateSurface(std::vector<int> idx, const DataView<double>& j) const;
+    double integrateSurface(std::vector<int> idx, const DataView<const double>& j) const;
 
-    double integrateVolume(int idx, double f) const;
-    double integrateVolume(int idx, const DataView<double>& f) const;
-    double integrateVolume(int idx, const DataView<const double>& f) const;
+    double integrateVolume(std::vector<int> idx, double f) const;
+    double integrateVolume(std::vector<int> idx, const DataView<double>& f) const;
+    double integrateVolume(std::vector<int> idx, const DataView<const double>& f) const;
 
     template<typename T>
     typename std::remove_const<T>::type interpolate(double x, const DataView<T>& f) const;
 
-    virtual double gradient(int idx, double f_lo, double f_hi) const = 0;
-    double gradient(int idx, const DataView<const double>& f) const;
-    double gradient(int idx, const DataView<double>& f) const;
+    virtual double gradient(std::vector<int> idx, double f_lo, double f_hi) const = 0;
+    double gradient(std::vector<int> idx, const DataView<const double>& f) const;
+    double gradient(std::vector<int> idx, const DataView<double>& f) const;
 
     bool operator==(const Mesh& other) const;
     bool operator!=(const Mesh& other) const;
 
     protected:
     double lower_;
-    int shape_; //!< Shape of the mesh
+    std::vector<int> shape_; //!< Shape of the mesh
     BoundaryType lower_bc_;
     BoundaryType upper_bc_;
     double step_; //!< Spacing between mesh points
-    int start_;
+    std::vector<int> start_;
 
     void validateBoundaryCondition() const;
 
diff --git a/src/mesh.cc b/src/mesh.cc
index 6e063f3..646f525 100644
--- a/src/mesh.cc
+++ b/src/mesh.cc
@@ -5,18 +5,28 @@
 namespace flyft
     {
 
-Mesh::Mesh(double lower_bound,
-           double upper_bound,
-           int shape,
+Mesh::Mesh(std::vector<double> lower_bound,
+           std::vector<double> upper_bound,
+           std::vector<int> shape,
            BoundaryType lower_bc,
            BoundaryType upper_bc)
     : lower_(lower_bound), shape_(shape), lower_bc_(lower_bc), upper_bc_(upper_bc), start_(0)
     {
-    step_ = (upper_bound - lower_bound) / shape_;
+    std::vector<double> step_;
+    std::transform(upper_bound.begin(),
+                   upper_bound.end(),
+                   lower_bound.begin(),
+                   step_.begin(),
+                   [](int i, int j) { return i - j });
+    std::transform(step_.begin(),
+                   step_.end(),
+                   shape_.begin(),
+                   step_.begin(),
+                   std::divides<double>());
     validateBoundaryCondition();
     }
 
-std::shared_ptr<Mesh> Mesh::slice(int start, int end) const
+std::shared_ptr<Mesh> Mesh::slice(std::vector<int> start, std::vector<int> end) const
     {
     if (lower_bc_ == BoundaryType::internal || upper_bc_ == BoundaryType::internal)
         {
@@ -37,14 +47,28 @@ std::shared_ptr<Mesh> Mesh::slice(int start, int end) const
     return m;
     }
 
-double Mesh::center(int i) const
+std::vector<double> Mesh::center(std::vector<int> i) const
     {
-    return lower_ + static_cast<double>(start_ + i + 0.5) * step_;
+    std::vector<double> temp;
+    std::transform(i.begin(),
+                   i.end(),
+                   start_.begin(),
+                   temp.begin(),
+                   [](int i, int j) { return static_cast<double>(i + j + 0.5) * step_ });
+    std::transform(temp.begin(), temp.end(), lower_.begin(), temp.begin(), std::plus<int>());
+    return temp;
     }
 
-int Mesh::bin(double x) const
+std::vector<int> Mesh::bin(std::vector<double> x) const
     {
-    return static_cast<int>((x - lower_) / step_) - start_;
+    std::vector<int> temp;
+    std::transform(x.begin(),
+                   x.end(),
+                   lower_.begin(),
+                   temp.begin(),
+                   [](int i, int j) { return static_cast<int>((i - j) / step_) });
+    std::transform(temp.begin(), temp.end(), start_.begin(), temp.begin().std::minus<int>);
+    return temp;
     }
 
 double Mesh::lower_bound() const
@@ -52,8 +76,15 @@ double Mesh::lower_bound() const
     return lower_bound(0);
     }
 
-double Mesh::lower_bound(int i) const
+std::vector<double> Mesh::lower_bound(std::vector<int> i) const
     {
+    std::vector<double> temp;
+    std::transform(i.begin(),
+                   i.end(),
+                   start_.begin(),
+                   temp.begin(),
+                   [](int i, int j) { return static_cast<double>((i + j) * step_) });
+    std::transform(temp.begin(), temp.end(), start_.begin(), temp.begin().std::minus<double>);
     return lower_ + (start_ + i) * step_;
     }
 
@@ -72,12 +103,12 @@ double Mesh::L() const
     return asLength(shape_);
     }
 
-int Mesh::shape() const
+std::vector<int> Mesh::shape() const
     {
     return shape_;
     }
 
-double Mesh::step() const
+std::vector<double> Mesh::step() const
     {
     return step_;
     }
@@ -92,34 +123,34 @@ double Mesh::asLength(int shape) const
     return shape * step_;
     }
 
-double Mesh::integrateSurface(int idx, double j_lo, double j_hi) const
+double Mesh::integrateSurface(std::vector<int> idx, double j_lo, double j_hi) const
     {
     return area(idx) * j_lo - area(idx + 1) * j_hi;
     }
 
-double Mesh::integrateSurface(int idx, const DataView<double>& j) const
+double Mesh::integrateSurface(std::vector<int> idx, const DataView<double>& j) const
     {
     return integrateSurface(idx, j(idx), j(idx + 1));
     }
 
-double Mesh::integrateSurface(int idx, const DataView<const double>& j) const
+double Mesh::integrateSurface(std::vector<int> idx, const DataView<const double>& j) const
     {
-    return integrateSurface(idx, j(idx), j(idx + 1));
+    return integrateSurface(std::vector<int>, j(idx), j(idx + 1));
     }
 
-double Mesh::integrateVolume(int idx, double f) const
+double Mesh::integrateVolume(std::vector<int> idx, double f) const
     {
     return volume(idx) * f;
     }
 
-double Mesh::integrateVolume(int idx, const DataView<double>& f) const
+double Mesh::integrateVolume(std::vector<int> idx, const DataView<double>& f) const
     {
-    return integrateVolume(idx, f(idx));
+    return integrateVolume(std::vector<int> idx, f(idx));
     }
 
-double Mesh::integrateVolume(int idx, const DataView<const double>& f) const
+double Mesh::integrateVolume(std::vector<int> idx, const DataView<const double>& f) const
     {
-    return integrateVolume(idx, f(idx));
+    return integrateVolume(std::vector<int> idx, f(idx));
     }
 
 double Mesh::gradient(int idx, const DataView<double>& f) const

From e3dabd943086eb298f21c7e413c108e7323a4c61 Mon Sep 17 00:00:00 2001
From: Mayukh Kundu <mzk0148@auburn.edu>
Date: Tue, 15 Apr 2025 13:19:00 -0500
Subject: [PATCH 2/6] Fix multi-dimensional indexing

---
 include/flyft/mesh.h | 55 ++++++++++++++++--------------
 src/mesh.cc          | 81 +++++++++++++++++++++++---------------------
 2 files changed, 72 insertions(+), 64 deletions(-)

diff --git a/include/flyft/mesh.h b/include/flyft/mesh.h
index ba70aff..56010f0 100644
--- a/include/flyft/mesh.h
+++ b/include/flyft/mesh.h
@@ -34,7 +34,7 @@ class Mesh
     std::vector<double> upper_bound() const;
 
     //! Get upper bound of bin
-    std::vector<double> upper_bound(int i) const;
+    std::vector<double> upper_bound(std::vector<int> i) const;
 
     //! Get surface area of lower edge of bin
     virtual double area(int i) const = 0;
@@ -49,7 +49,7 @@ class Mesh
     std::vector<int> bin(double x) const;
 
     //! Length of the mesh
-    double L() const;
+    std::vector<double> L() const;
 
     //! Shape of the mesh
     std::vector<int> shape() const;
@@ -63,9 +63,9 @@ class Mesh
     //! Boundary condition on upper bound of mesh
     BoundaryType upper_boundary_condition() const;
 
-    int asShape(double dx) const;
+    std::vector<int> asShape(std::vector<double> dx) const;
 
-    double asLength(int shape) const;
+    std::vector<double> asLength(std::vector<int> shape) const;
 
     double integrateSurface(std::vector<int> idx, double j_lo, double j_hi) const;
     double integrateSurface(std::vector<int> idx, const DataView<double>& j) const;
@@ -86,11 +86,11 @@ class Mesh
     bool operator!=(const Mesh& other) const;
 
     protected:
-    double lower_;
+    std::vector<double> lower_;
     std::vector<int> shape_; //!< Shape of the mesh
     BoundaryType lower_bc_;
     BoundaryType upper_bc_;
-    double step_; //!< Spacing between mesh points
+    std::vector<double> step_; //!< Spacing between mesh points
     std::vector<int> start_;
 
     void validateBoundaryCondition() const;
@@ -99,28 +99,33 @@ class Mesh
     };
 
 template<typename T>
-typename std::remove_const<T>::type Mesh::interpolate(double x, const DataView<T>& f) const
+typename std::remove_const<T>::type Mesh::interpolate(std::vector<double> x,
+                                                      const DataView<T>& f) const
     {
-    const auto idx = bin(x);
-    const auto x_c = center(idx);
-    double x_0, x_1;
-    typename std::remove_const<T>::type f_0, f_1;
-    if (x < x_c)
+    std::vector<T> temp;
+    for (int k = 0; k < shape_.size(); ++k)
         {
-        x_0 = center(idx - 1);
-        x_1 = x_c;
-        f_0 = f(idx - 1);
-        f_1 = f(idx);
+        const auto idx = bin(x[k]);
+        const auto x_c = center(idx);
+        double x_0, x_1;
+        typename std::remove_const<T>::type f_0, f_1;
+        if (x < x_c)
+            {
+            x_0 = center(idx - 1);
+            x_1 = x_c;
+            f_0 = f(idx - 1);
+            f_1 = f(idx);
+            }
+        else
+            {
+            x_0 = x_c;
+            x_1 = center(idx + 1);
+            f_0 = f(idx);
+            f_1 = f(idx + 1);
+            }
+        temp.push_back(f_0 + (x - x_0) * (f_1 - f_0) / (x_1 - x_0));
         }
-    else
-        {
-        x_0 = x_c;
-        x_1 = center(idx + 1);
-        f_0 = f(idx);
-        f_1 = f(idx + 1);
-        }
-
-    return f_0 + (x - x_0) * (f_1 - f_0) / (x_1 - x_0);
+    return temp;
     }
 
     } // namespace flyft
diff --git a/src/mesh.cc b/src/mesh.cc
index 646f525..1a081fe 100644
--- a/src/mesh.cc
+++ b/src/mesh.cc
@@ -13,16 +13,10 @@ Mesh::Mesh(std::vector<double> lower_bound,
     : lower_(lower_bound), shape_(shape), lower_bc_(lower_bc), upper_bc_(upper_bc), start_(0)
     {
     std::vector<double> step_;
-    std::transform(upper_bound.begin(),
-                   upper_bound.end(),
-                   lower_bound.begin(),
-                   step_.begin(),
-                   [](int i, int j) { return i - j });
-    std::transform(step_.begin(),
-                   step_.end(),
-                   shape_.begin(),
-                   step_.begin(),
-                   std::divides<double>());
+    for (int idx = 0; idx < shape_.size(); ++idx)
+        {
+        step_.push_back((upper_bound[idx] - lower_bound[idx]) / shape_[idx]);
+        }
     validateBoundaryCondition();
     }
 
@@ -50,55 +44,54 @@ std::shared_ptr<Mesh> Mesh::slice(std::vector<int> start, std::vector<int> end)
 std::vector<double> Mesh::center(std::vector<int> i) const
     {
     std::vector<double> temp;
-    std::transform(i.begin(),
-                   i.end(),
-                   start_.begin(),
-                   temp.begin(),
-                   [](int i, int j) { return static_cast<double>(i + j + 0.5) * step_ });
-    std::transform(temp.begin(), temp.end(), lower_.begin(), temp.begin(), std::plus<int>());
+    for (int idx; idx < shape_.size(); ++idx)
+        {
+        temp.push_back(lower_[idx] + static_cast<double>(start_[idx] + i[idx] + 0.5))
+        }
     return temp;
     }
 
 std::vector<int> Mesh::bin(std::vector<double> x) const
     {
     std::vector<int> temp;
-    std::transform(x.begin(),
-                   x.end(),
-                   lower_.begin(),
-                   temp.begin(),
-                   [](int i, int j) { return static_cast<int>((i - j) / step_) });
-    std::transform(temp.begin(), temp.end(), start_.begin(), temp.begin().std::minus<int>);
+    for (int idx; idx < shape_.size(); ++idx)
+        {
+        temp.push_back(static_cast<int>((x[idx] - lower_[idx]) / step_[idx]) - start_[idx])
+        }
     return temp;
     }
 
-double Mesh::lower_bound() const
+std::vector<double> Mesh::lower_bound() const
     {
-    return lower_bound(0);
+    std::vector<int> temp(shape_.size(), 0) return lower_bound(temp);
     }
 
 std::vector<double> Mesh::lower_bound(std::vector<int> i) const
     {
     std::vector<double> temp;
-    std::transform(i.begin(),
-                   i.end(),
-                   start_.begin(),
-                   temp.begin(),
-                   [](int i, int j) { return static_cast<double>((i + j) * step_) });
-    std::transform(temp.begin(), temp.end(), start_.begin(), temp.begin().std::minus<double>);
-    return lower_ + (start_ + i) * step_;
+    for (int idx; idx < shape_.size(); ++idx)
+        {
+        temp.push_back(lower_[idx] + (start_[idx] + i[idx]) * step_[idx]);
+        }
+    return temp
     }
 
-double Mesh::upper_bound() const
+std::vector<double> Mesh::upper_bound() const
     {
     return lower_bound(shape_);
     }
 
-double Mesh::upper_bound(int i) const
+std::vector<double> Mesh::upper_bound(std::vector<int> i) const
     {
-    return lower_bound(i + 1);
+    std::vector<double> temp;
+    for (int idx; idx < shape_.size(); ++idx)
+        {
+        temp.push_back(i[idx] + 1);
+        }
+    return lower_bound(temp);
     }
 
-double Mesh::L() const
+std::vector<double> Mesh::L() const
     {
     return asLength(shape_);
     }
@@ -113,14 +106,24 @@ std::vector<double> Mesh::step() const
     return step_;
     }
 
-int Mesh::asShape(double dx) const
+std::vector<int> Mesh::asShape(std::vector<double> dx) const
     {
-    return static_cast<int>(std::ceil(dx / step_));
+    std::vector<int> temp;
+    for (int idx = 0; idx < shape_.size(); ++idx)
+        {
+        temp.push_back(static_cast<int>(std::ceil(dx[idx] / step_[idx])));
+        }
+    return temp;
     }
 
-double Mesh::asLength(int shape) const
+std::vector<double> Mesh::asLength(std::vector<int> shape) const
     {
-    return shape * step_;
+    std::vector<int> temp;
+    for (int idx = 0; idx < shape_.size(); ++idx)
+        {
+        temp.push_back(shape[idx] * step_[idx]);
+        }
+    return temp;
     }
 
 double Mesh::integrateSurface(std::vector<int> idx, double j_lo, double j_hi) const

From f374be9d6ef27ca67072fd0391f322c6edb70cf7 Mon Sep 17 00:00:00 2001
From: Mayukh Kundu <mzk0148@auburn.edu>
Date: Wed, 16 Apr 2025 16:17:43 -0500
Subject: [PATCH 3/6] Apply fix and add documentation

---
 include/flyft/mesh.h | 132 +++++++++++++++++++++++++++++++++++--------
 src/mesh.cc          |  51 +++++++++--------
 2 files changed, 136 insertions(+), 47 deletions(-)

diff --git a/include/flyft/mesh.h b/include/flyft/mesh.h
index 56010f0..f41929f 100644
--- a/include/flyft/mesh.h
+++ b/include/flyft/mesh.h
@@ -8,47 +8,104 @@
 
 namespace flyft
     {
-
+//! Multidimensional mesh
+/*!
+ * Multidimensional mesh maps N-dimensional indexing for
+ * the anistropic particles into one-dimensional flux.
+ */
 class Mesh
     {
     public:
+    //! No default constructor
     Mesh() = delete;
+
+    //! Constructor
+    /*!
+     * \param lower_bound Lower bound of the N-dimensional mesh.
+     * \param upper_bound Upper bound of the N-dimensional mesh.
+     * \param shape Shape of the index array.
+     * \param lower_bc Boundary condition of the lower bound.
+     * \param upper_bc Boundary condition of the upper bound.
+     */
     Mesh(double lower_bound,
          double upper_bound,
          int shape,
          BoundaryType lower_bc,
          BoundaryType upper_bc);
 
-    std::shared_ptr<Mesh> slice(std::vector<int> start, std::vector<int> end) const;
+    //! Slice multidimensional index
+    /*!
+     * \param start Multidimensional start index.
+     * \param end Multidimensional end index.
+     * \return One-dimensional index.
+     *
+     */
+    std::shared_ptr<Mesh> slice(const std::vector<int>& start, const std::vector<int>& end) const;
 
     //! Get position on the mesh, defined as center of bin
-    std::vector<double> center(std::vector<int> i) const;
+    /*!
+     * \param i Multidimensional bin index
+     * \return Center of the multidimensional bin
+     */
+    std::vector<double> center(const std::vector<int>& i) const;
 
     //! Lower bound of entire mesh
+    /*!
+     * \return Lower bound of the multidimensional mesh
+     */
     std::vector<double> lower_bound() const;
 
     //! Get lower bound of bin
-    std::vector<double> lower_bound(std::vector<int> i) const;
+    /*!
+     * \param i Multidimensional bin index
+     * \return Lower bound of the multidimensional bin
+     */
+    std::vector<double> lower_bound(const std::vector<int>& i) const;
 
     //! Upper bound of entire mesh
+    /*!
+     * \return Upper bound of the multidimensional mesh
+     */
     std::vector<double> upper_bound() const;
 
     //! Get upper bound of bin
-    std::vector<double> upper_bound(std::vector<int> i) const;
+    /*!
+     * \param i Multidimensional bin index
+     * \return Upper bound of the multidimensional bin
+     */
+    std::vector<double> upper_bound(const std::vector<int>& i) const;
 
     //! Get surface area of lower edge of bin
+    /*!
+     * \param i Multidimensional bin index
+     * \return Upper bound of the multidimensional bin
+     */
     virtual double area(int i) const = 0;
 
     //! Get volume of mesh
+    /*!
+     * \return Volume of the mesh
+     */
     virtual double volume() const = 0;
 
     //! Get volume of bin
+    /*!
+     * \param i Multidimensional bin index
+     * \return Volume of the bin
+     */
     virtual double volume(int i) const = 0;
 
     //! Get the bin for a coordinate
+    /*!
+     * \param x Coordinate position of the bin
+     * \return Bin index
+     */
     std::vector<int> bin(double x) const;
 
     //! Length of the mesh
+    /*!
+     * \return Length of the mesh
+     */
     std::vector<double> L() const;
 
     //! Shape of the mesh
@@ -63,33 +120,62 @@ class Mesh
     //! Boundary condition on upper bound of mesh
     BoundaryType upper_boundary_condition() const;
 
-    std::vector<int> asShape(std::vector<double> dx) const;
-
-    std::vector<double> asLength(std::vector<int> shape) const;
-
-    double integrateSurface(std::vector<int> idx, double j_lo, double j_hi) const;
-    double integrateSurface(std::vector<int> idx, const DataView<double>& j) const;
-    double integrateSurface(std::vector<int> idx, const DataView<const double>& j) const;
-
-    double integrateVolume(std::vector<int> idx, double f) const;
-    double integrateVolume(std::vector<int> idx, const DataView<double>& f) const;
-    double integrateVolume(std::vector<int> idx, const DataView<const double>& f) const;
-
+    //! Get the start index of the mesh
+    /*!
+     * \param dx Step size of the mesh
+     * \return Start index of the mesh
+     */
+    std::vector<int> asShape(const std::vector<double>& dx) const;
+
+    //! Get the length of the mesh
+    /*!
+     * \param shape Shape of the mesh
+     * \return Length of the mesh
+     */
+    std::vector<double> asLength(const std::vector<int>& shape) const;
+
+    //! Get the integral of the surface over the mesh
+    /*!
+     * \param shape Shape of the mesh
+     * \return Length of the mesh
+     */
+    double integrateSurface(const std::vector<int>& idx, double j_lo, double j_hi) const;
+    double integrateSurface(const std::vector<int>& idx, const DataView<double>& j) const;
+    double integrateSurface(const std::vector<int>& idx, const DataView<const double>& j) const;
+    //! Get the integral of the volume over the mesh
+    /*!
+     * \param idx Multidimensional index
+     * \param f Function to integrate
+     * \return Integral of the function over the mesh
+     */
+    double integrateVolume(const std::vector<int>& idx, double f) const;
+    double integrateVolume(const std::vector<int>& idx, const DataView<double>& f) const;
+    double integrateVolume(const std::vector<int>& idx, const DataView<const double>& f) const;
+
+    //! Interpolate a function on the mesh
+    /*!
+     * \param x Coordinate position of the bin
+     * \param f Function to interpolate
+     * \return Interpolated function value
+     *
+     * The function is interpolated using linear interpolation.
+     */
     template<typename T>
     typename std::remove_const<T>::type interpolate(double x, const DataView<T>& f) const;
 
-    virtual double gradient(std::vector<int> idx, double f_lo, double f_hi) const = 0;
-    double gradient(std::vector<int> idx, const DataView<const double>& f) const;
-    double gradient(std::vector<int> idx, const DataView<double>& f) const;
+    //! Get the gradient of the mesh
+    virtual double gradient(const std::vector<int>& idx, double f_lo, double f_hi) const = 0;
+    double gradient(const std::vector<int>& idx, const DataView<const double>& f) const;
+    double gradient(const std::vector<int>& idx, const DataView<double>& f) const;
 
     bool operator==(const Mesh& other) const;
     bool operator!=(const Mesh& other) const;
 
     protected:
     std::vector<double> lower_;
-    std::vector<int> shape_; //!< Shape of the mesh
-    BoundaryType lower_bc_;
-    BoundaryType upper_bc_;
+    std::vector<int> shape_;   //!< Shape of the mesh
+    BoundaryType lower_bc_;    //!< Boundary condition of the lower bound
+    BoundaryType upper_bc_;    //!< Boundary condition of the upper bound
     std::vector<double> step_; //!< Spacing between mesh points
     std::vector<int> start_;
 
diff --git a/src/mesh.cc b/src/mesh.cc
index 1a081fe..1f66562 100644
--- a/src/mesh.cc
+++ b/src/mesh.cc
@@ -4,7 +4,6 @@
 
 namespace flyft
     {
-
 Mesh::Mesh(std::vector<double> lower_bound,
            std::vector<double> upper_bound,
            std::vector<int> shape,
@@ -20,7 +19,7 @@ Mesh::Mesh(std::vector<double> lower_bound,
     validateBoundaryCondition();
     }
 
-std::shared_ptr<Mesh> Mesh::slice(std::vector<int> start, std::vector<int> end) const
+std::shared_ptr<Mesh> Mesh::slice(const std::vector<int>& start, const std::vector<int>& end) const
     {
     if (lower_bc_ == BoundaryType::internal || upper_bc_ == BoundaryType::internal)
         {
@@ -28,20 +27,24 @@ std::shared_ptr<Mesh> Mesh::slice(std::vector<int> start, std::vector<int> end)
         }
 
     auto m = clone();
-    m->start_ = start;
-    m->shape_ = end - start;
-    if (start > 0)
-        {
-        m->lower_bc_ = BoundaryType::internal;
-        }
-    if (end < shape_)
+    for (int idx = 0; idx < start.size(); ++idx)
         {
-        m->upper_bc_ = BoundaryType::internal;
+        m->start_[idx] = start[idx];
+        m->shape_[idx] = end[idx] - start[idx];
+        if (start[idx] > 0)
+            {
+            m->lower_bc_ = BoundaryType::internal;
+            }
+        if (end < shape_)
+            {
+            m->upper_bc_ = BoundaryType::internal;
+            }
         }
+
     return m;
     }
 
-std::vector<double> Mesh::center(std::vector<int> i) const
+std::vector<double> Mesh::center(const std::vector<int>& i) const
     {
     std::vector<double> temp;
     for (int idx; idx < shape_.size(); ++idx)
@@ -51,7 +54,7 @@ std::vector<double> Mesh::center(std::vector<int> i) const
     return temp;
     }
 
-std::vector<int> Mesh::bin(std::vector<double> x) const
+std::vector<int> Mesh::bin(const std::vector<double>& x) const
     {
     std::vector<int> temp;
     for (int idx; idx < shape_.size(); ++idx)
@@ -66,7 +69,7 @@ std::vector<double> Mesh::lower_bound() const
     std::vector<int> temp(shape_.size(), 0) return lower_bound(temp);
     }
 
-std::vector<double> Mesh::lower_bound(std::vector<int> i) const
+std::vector<double> Mesh::lower_bound(const std::vector<int>& i) const
     {
     std::vector<double> temp;
     for (int idx; idx < shape_.size(); ++idx)
@@ -81,7 +84,7 @@ std::vector<double> Mesh::upper_bound() const
     return lower_bound(shape_);
     }
 
-std::vector<double> Mesh::upper_bound(std::vector<int> i) const
+std::vector<double> Mesh::upper_bound(const std::vector<int>& i) const
     {
     std::vector<double> temp;
     for (int idx; idx < shape_.size(); ++idx)
@@ -106,7 +109,7 @@ std::vector<double> Mesh::step() const
     return step_;
     }
 
-std::vector<int> Mesh::asShape(std::vector<double> dx) const
+std::vector<int> Mesh::asShape(const std::vector<double>& dx) const
     {
     std::vector<int> temp;
     for (int idx = 0; idx < shape_.size(); ++idx)
@@ -116,7 +119,7 @@ std::vector<int> Mesh::asShape(std::vector<double> dx) const
     return temp;
     }
 
-std::vector<double> Mesh::asLength(std::vector<int> shape) const
+std::vector<double> Mesh::asLength(const std::vector<int>& shape) const
     {
     std::vector<int> temp;
     for (int idx = 0; idx < shape_.size(); ++idx)
@@ -126,19 +129,19 @@ std::vector<double> Mesh::asLength(std::vector<int> shape) const
     return temp;
     }
 
-double Mesh::integrateSurface(std::vector<int> idx, double j_lo, double j_hi) const
+double Mesh::integrateSurface(const std::vector<int>& idx, double j_lo, double j_hi) const
     {
     return area(idx) * j_lo - area(idx + 1) * j_hi;
     }
 
-double Mesh::integrateSurface(std::vector<int> idx, const DataView<double>& j) const
+double Mesh::integrateSurface(const std::vector<int>& idx, const DataView<double>& j) const
     {
     return integrateSurface(idx, j(idx), j(idx + 1));
     }
 
-double Mesh::integrateSurface(std::vector<int> idx, const DataView<const double>& j) const
+double Mesh::integrateSurface(const std::vector<int>& idx, const DataView<const double>& j) const
     {
-    return integrateSurface(std::vector<int>, j(idx), j(idx + 1));
+    return integrateSurface(const std::vector<int>& idx, j(idx), j(idx + 1));
     }
 
 double Mesh::integrateVolume(std::vector<int> idx, double f) const
@@ -146,14 +149,14 @@ double Mesh::integrateVolume(std::vector<int> idx, double f) const
     return volume(idx) * f;
     }
 
-double Mesh::integrateVolume(std::vector<int> idx, const DataView<double>& f) const
+double Mesh::integrateVolume(const std::vector<int>& idx, const DataView<double>& f) const
     {
-    return integrateVolume(std::vector<int> idx, f(idx));
+    return integrateVolume(const std::vector<int>& idx, f(idx));
     }
 
-double Mesh::integrateVolume(std::vector<int> idx, const DataView<const double>& f) const
+double Mesh::integrateVolume(const std::vector<int>& idx, const DataView<const double>& f) const
     {
-    return integrateVolume(std::vector<int> idx, f(idx));
+    return integrateVolume(const std::vector<int>& idx, f(idx));
     }
 
 double Mesh::gradient(int idx, const DataView<double>& f) const

From 1ece6abfd4c288c3e0e9fef85e9353d6980c391a Mon Sep 17 00:00:00 2001
From: Mayukh Kundu <mzk0148@auburn.edu>
Date: Thu, 17 Apr 2025 22:23:58 -0500
Subject: [PATCH 4/6] Minor fix

---
 src/mesh.cc | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/mesh.cc b/src/mesh.cc
index 1f66562..deae9b3 100644
--- a/src/mesh.cc
+++ b/src/mesh.cc
@@ -11,7 +11,6 @@ Mesh::Mesh(std::vector<double> lower_bound,
            BoundaryType upper_bc)
     : lower_(lower_bound), shape_(shape), lower_bc_(lower_bc), upper_bc_(upper_bc), start_(0)
     {
-    std::vector<double> step_;
     for (int idx = 0; idx < shape_.size(); ++idx)
         {
         step_.push_back((upper_bound[idx] - lower_bound[idx]) / shape_[idx]);

From fc8090e68ba3c9ed760651936a1732cdbd20215c Mon Sep 17 00:00:00 2001
From: Mayukh Kundu <mzk0148@auburn.edu>
Date: Thu, 5 Jun 2025 11:25:17 -0500
Subject: [PATCH 5/6] Apply review suggestions and minor fixes

---
 .pre-commit-config.yaml |  2 +-
 include/flyft/mesh.h    | 48 +++++++++++++++++++----------------------
 src/mesh.cc             |  7 +++---
 3 files changed, 27 insertions(+), 30 deletions(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 1b63521..9f0779d 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -21,7 +21,7 @@ repos:
     hooks:
       - id: isort
   - repo: https://github.com/psf/black-pre-commit-mirror
-    rev: 24.10.0
+    rev: 24.8.0
     hooks:
       - id: black
   - repo: https://github.com/PyCQA/flake8
diff --git a/include/flyft/mesh.h b/include/flyft/mesh.h
index f41929f..6a24370 100644
--- a/include/flyft/mesh.h
+++ b/include/flyft/mesh.h
@@ -93,14 +93,14 @@ class Mesh
      * \param i Multidimensional bin index
      * \return Volume of the bin
      */
-    virtual double volume(int i) const = 0;
+    virtual std::vector<double> volume(int i) const = 0;
 
     //! Get the bin for a coordinate
     /*!
      * \param x Coordinate position of the bin
      * \return Bin index
      */
-    std::vector<int> bin(double x) const;
+    std::vector<int> bin(const std::vector<double>& x) const;
 
     //! Length of the mesh
     /*!
@@ -142,6 +142,7 @@ class Mesh
     double integrateSurface(const std::vector<int>& idx, double j_lo, double j_hi) const;
     double integrateSurface(const std::vector<int>& idx, const DataView<double>& j) const;
     double integrateSurface(const std::vector<int>& idx, const DataView<const double>& j) const;
+
     //! Get the integral of the volume over the mesh
     /*!
      * \param idx Multidimensional index
@@ -185,33 +186,28 @@ class Mesh
     };
 
 template<typename T>
-typename std::remove_const<T>::type Mesh::interpolate(std::vector<double> x,
-                                                      const DataView<T>& f) const
+typename std::remove_const<T>::type Mesh::interpolate(double x, const DataView<T>& f) const
     {
-    std::vector<T> temp;
-    for (int k = 0; k < shape_.size(); ++k)
+    const auto idx = bin(x);
+    const auto x_c = center(idx);
+    double x_0, x_1;
+    typename std::remove_const<T>::type f_0, f_1;
+    if (x < x_c)
+        {
+        x_0 = center(idx - 1);
+        x_1 = x_c;
+        f_0 = f(idx - 1);
+        f_1 = f(idx);
+        }
+    else
         {
-        const auto idx = bin(x[k]);
-        const auto x_c = center(idx);
-        double x_0, x_1;
-        typename std::remove_const<T>::type f_0, f_1;
-        if (x < x_c)
-            {
-            x_0 = center(idx - 1);
-            x_1 = x_c;
-            f_0 = f(idx - 1);
-            f_1 = f(idx);
-            }
-        else
-            {
-            x_0 = x_c;
-            x_1 = center(idx + 1);
-            f_0 = f(idx);
-            f_1 = f(idx + 1);
-            }
-        temp.push_back(f_0 + (x - x_0) * (f_1 - f_0) / (x_1 - x_0));
+        x_0 = x_c;
+        x_1 = center(idx + 1);
+        f_0 = f(idx);
+        f_1 = f(idx + 1);
         }
-    return temp;
+
+    return f_0 + (x - x_0) * (f_1 - f_0) / (x_1 - x_0);
     }
 
     } // namespace flyft
diff --git a/src/mesh.cc b/src/mesh.cc
index deae9b3..3303edc 100644
--- a/src/mesh.cc
+++ b/src/mesh.cc
@@ -48,7 +48,7 @@ std::vector<double> Mesh::center(const std::vector<int>& i) const
     std::vector<double> temp;
     for (int idx; idx < shape_.size(); ++idx)
         {
-        temp.push_back(lower_[idx] + static_cast<double>(start_[idx] + i[idx] + 0.5))
+        temp.push_back(lower_[idx] + static_cast<double>(start_[idx] + i[idx] + 0.5));
         }
     return temp;
     }
@@ -58,14 +58,15 @@ std::vector<int> Mesh::bin(const std::vector<double>& x) const
     std::vector<int> temp;
     for (int idx; idx < shape_.size(); ++idx)
         {
-        temp.push_back(static_cast<int>((x[idx] - lower_[idx]) / step_[idx]) - start_[idx])
+        temp.push_back(static_cast<int>((x[idx] - lower_[idx]) / step_[idx]) - start_[idx]);
         }
     return temp;
     }
 
 std::vector<double> Mesh::lower_bound() const
     {
-    std::vector<int> temp(shape_.size(), 0) return lower_bound(temp);
+    std::vector<int> temp(shape_.size(), 0);
+    return lower_bound(temp);
     }
 
 std::vector<double> Mesh::lower_bound(const std::vector<int>& i) const

From 361bd9c1e3aa75bcd5aefc4aa6388a3266ba46d0 Mon Sep 17 00:00:00 2001
From: Mayukh Kundu <79770572+mayukh33@users.noreply.github.com>
Date: Thu, 5 Jun 2025 11:47:55 -0500
Subject: [PATCH 6/6] Revert .pre-commit-config.yaml

---
 .pre-commit-config.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 9f0779d..1b63521 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -21,7 +21,7 @@ repos:
     hooks:
       - id: isort
   - repo: https://github.com/psf/black-pre-commit-mirror
-    rev: 24.8.0
+    rev: 24.10.0
     hooks:
       - id: black
   - repo: https://github.com/PyCQA/flake8