diff --git a/include/flyft/mesh.h b/include/flyft/mesh.h
index ab85cc0..6a24370 100644
--- a/include/flyft/mesh.h
+++ b/include/flyft/mesh.h
@@ -8,54 +8,111 @@
 
 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(int start, 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
-    double center(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
-    double lower_bound() const;
+    /*!
+     * \return Lower bound of the multidimensional mesh
+     */
+    std::vector<double> lower_bound() const;
 
     //! Get lower bound of bin
-    double lower_bound(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
-    double upper_bound() const;
+    /*!
+     * \return Upper bound of the multidimensional mesh
+     */
+    std::vector<double> upper_bound() const;
 
     //! Get upper bound of bin
-    double upper_bound(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
-    virtual double volume(int i) const = 0;
+    /*!
+     * \param i Multidimensional bin index
+     * \return Volume of the bin
+     */
+    virtual std::vector<double> volume(int i) const = 0;
 
     //! Get the bin for a coordinate
-    int bin(double x) const;
+    /*!
+     * \param x Coordinate position of the bin
+     * \return Bin index
+     */
+    std::vector<int> bin(const std::vector<double>& x) const;
 
     //! Length of the mesh
-    double L() const;
+    /*!
+     * \return Length of the mesh
+     */
+    std::vector<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;
@@ -63,35 +120,65 @@ class Mesh
     //! Boundary condition on upper bound of mesh
     BoundaryType upper_boundary_condition() const;
 
-    int asShape(double dx) const;
-
-    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 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;
-
+    //! 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(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;
+    //! 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:
-    double lower_;
-    int shape_; //!< Shape of the mesh
-    BoundaryType lower_bc_;
-    BoundaryType upper_bc_;
-    double step_; //!< Spacing between mesh points
-    int start_;
+    std::vector<double> lower_;
+    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_;
 
     void validateBoundaryCondition() const;
 
diff --git a/src/mesh.cc b/src/mesh.cc
index 6e063f3..3303edc 100644
--- a/src/mesh.cc
+++ b/src/mesh.cc
@@ -4,19 +4,21 @@
 
 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_;
+    for (int idx = 0; idx < shape_.size(); ++idx)
+        {
+        step_.push_back((upper_bound[idx] - lower_bound[idx]) / shape_[idx]);
+        }
     validateBoundaryCondition();
     }
 
-std::shared_ptr<Mesh> Mesh::slice(int start, 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)
         {
@@ -24,102 +26,137 @@ std::shared_ptr<Mesh> Mesh::slice(int start, int end) const
         }
 
     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;
     }
 
-double Mesh::center(int i) const
+std::vector<double> Mesh::center(const std::vector<int>& i) const
     {
-    return lower_ + static_cast<double>(start_ + i + 0.5) * step_;
+    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));
+        }
+    return temp;
     }
 
-int Mesh::bin(double x) const
+std::vector<int> Mesh::bin(const std::vector<double>& x) const
     {
-    return static_cast<int>((x - lower_) / step_) - start_;
+    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]);
+        }
+    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);
     }
 
-double Mesh::lower_bound(int i) const
+std::vector<double> Mesh::lower_bound(const std::vector<int>& i) const
     {
-    return lower_ + (start_ + i) * step_;
+    std::vector<double> temp;
+    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(const 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_);
     }
 
-int Mesh::shape() const
+std::vector<int> Mesh::shape() const
     {
     return shape_;
     }
 
-double Mesh::step() const
+std::vector<double> Mesh::step() const
     {
     return step_;
     }
 
-int Mesh::asShape(double dx) const
+std::vector<int> Mesh::asShape(const 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(const 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(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(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(int idx, const DataView<const double>& j) const
+double Mesh::integrateSurface(const std::vector<int>& idx, const DataView<const double>& j) const
     {
-    return integrateSurface(idx, j(idx), j(idx + 1));
+    return integrateSurface(const std::vector<int>& idx, 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(const std::vector<int>& idx, const DataView<double>& f) const
     {
-    return integrateVolume(idx, f(idx));
+    return integrateVolume(const std::vector<int>& idx, f(idx));
     }
 
-double Mesh::integrateVolume(int idx, const DataView<const double>& f) const
+double Mesh::integrateVolume(const std::vector<int>& idx, const DataView<const double>& f) const
     {
-    return integrateVolume(idx, f(idx));
+    return integrateVolume(const std::vector<int>& idx, f(idx));
     }
 
 double Mesh::gradient(int idx, const DataView<double>& f) const