diff --git a/include/PolarGrid/multiindex.h b/include/PolarGrid/multiindex.h
deleted file mode 100755
index 32508f0b..00000000
--- a/include/PolarGrid/multiindex.h
+++ /dev/null
@@ -1,46 +0,0 @@
-#pragma once
-
-#include <algorithm>
-#include <cassert>
-#include <iterator>
-
-#include "../common/space_dimension.h"
-
-class MultiIndex
-{
-public:
-    static_assert(space_dimension > 0 && space_dimension <= 3, "Invalid space dimension");
-
-    //! Creates a multi index with undefined coordinate values.
-    MultiIndex() = default;
-
-    //! Initializes all coordinate values of the multi index with `value`.
-    explicit MultiIndex(int value);
-    //! Initializes 2-dimensional multi index coordinate values with `i` and `j`.
-    //! Works only if `space_dimension` is 2.
-    explicit MultiIndex(int i, int j);
-    //! Initializes 3-dimensional multi index coordinate values with `i`, `j` and `k`.
-    //! Works only if `space_dimension` is 3.
-    explicit MultiIndex(int i, int j, int k);
-
-    //! Returns the size of the multi index.
-    //! This is equal to `space_dimension`.
-    int size() const;
-
-    //! Tests if two multi indices are equal
-    bool operator==(const MultiIndex& other) const;
-    //! Tests if two multi indices are inequal
-    bool operator!=(const MultiIndex& other) const;
-
-    //! Returns the `i`th coordinate value of the multi index.
-    const int& operator[](int i) const;
-
-    //! Returns a mutable reference to the `i`th coordinate value of the multi index.
-    int& operator[](int i);
-
-    int* data();
-    const int* data() const;
-
-private:
-    int data_[space_dimension];
-};
\ No newline at end of file
diff --git a/include/PolarGrid/point.h b/include/PolarGrid/point.h
deleted file mode 100755
index 2a288e9c..00000000
--- a/include/PolarGrid/point.h
+++ /dev/null
@@ -1,45 +0,0 @@
-#pragma once
-
-#include <algorithm>
-#include <cassert>
-#include <iterator>
-
-#include "../common/equals.h"
-#include "../common/space_dimension.h"
-
-class Point
-{
-public:
-    static_assert(space_dimension > 0 && space_dimension <= 3, "Invalid space dimension");
-
-    //! Creates a point with undefined coordinate values.
-    Point() = default;
-
-    //! Initializes all coordinate values of the point with `value`.
-    explicit Point(double value);
-    //! Initializes 2-dimensional point coordinate values with `x` and `y`.
-    //! Works only if `space_dimension` is 2.
-    explicit Point(double x, double y);
-    //! Initializes 3-dimensional point coordinate values with `x`, `y` and `z`.
-    //! Works only if `space_dimension` is 3.
-    explicit Point(double x, double y, double z);
-
-    //! Returns the size of the point.
-    //! This is equal to `space_dimension`.
-    int size() const;
-
-    //! Returns the `i`th coordinate value of the point.
-    double operator[](int i) const;
-
-    //! Returns a mutable reference to the `i`th coordinate value of the point.
-    double& operator[](int i);
-
-private:
-    double data_[space_dimension];
-};
-
-bool equals(const Point& lhs, const Point& rhs);
-double norm(const Point& point);
-void add(Point& result, const Point& lhs, const Point& rhs);
-void subtract(Point& result, const Point& lhs, const Point& rhs);
-void multiply(Point& result, double scalar, const Point& lhs);
\ No newline at end of file
diff --git a/include/PolarGrid/polargrid.h b/include/PolarGrid/polargrid.h
index 15b9e05d..a05c44d9 100644
--- a/include/PolarGrid/polargrid.h
+++ b/include/PolarGrid/polargrid.h
@@ -19,13 +19,6 @@
 
 #include "../common/equals.h"
 
-#include "../PolarGrid/multiindex.h"
-#include "../PolarGrid/point.h"
-
-// The PolarGrid class implements a donut-shaped 2D grid.
-// It is designed to handle polar coordinates, providing functionalities
-// for storing data points and performing operations on them.
-
 class PolarGrid
 {
 public:
@@ -44,12 +37,11 @@ class PolarGrid
     // Optimized, inlined indexing.
     int wrapThetaIndex(const int unwrapped_theta_index) const;
     int index(const int r_index, const int unwrapped_theta_index) const;
-    int fastIndex(const int r_index, const int theta_index) const;
     void multiIndex(const int node_index, int& r_index, int& theta_index) const;
 
+    // Grid Parameters
     // Number of grid nodes
     int numberOfNodes() const;
-    // Grid Parameters
     // Get the number of grid points in radial direction
     int nr() const;
     // Get the number of angular divisions
@@ -58,9 +50,6 @@ class PolarGrid
     double radius(const int r_index) const;
     // Get the angle at a specific angular index
     double theta(const int theta_index) const;
-    // Get all radii and angles available which define the grid
-    const std::vector<double>& radii() const;
-    const std::vector<double>& angles() const;
 
     // Grid distances
     // Get the radial distance to the next consecutive radial node at a specified radial index.
@@ -80,48 +69,8 @@ class PolarGrid
     // Get the number of nodes in radial smoother.
     int numberRadialSmootherNodes() const;
 
-    // ------------------------------------------- //
-    // Unoptimized Indexing and Neighbor Retrieval //
-    // ------------------------------------------- //
-    // Node Indexing (based on the combined circle-radial smoother)
-    // Get the index of a node based on its position.
-    int index(const MultiIndex& position) const;
-    // Get the position of a node based on its index.
-    MultiIndex multiIndex(const int node_index) const;
-    // Get the polar coordinates of a node based on its position.
-    Point polarCoordinates(const MultiIndex& position) const;
-    // Get adjacent neighbors of a node.
-    // If the neighbor index is -1, then there is no neighboring node in that direction.
-    // - The first entry (neighbors[0]) represents the radial direction (r):
-    //   - First: inward neighbor (r - 1)
-    //   - Second: outward neighbor (r + 1)
-    // - The second entry (neighbors[1]) represents the angular direction (theta):
-    //   - First: counterclockwise neighbor (theta - 1)
-    //   - Second: clockwise neighbor (theta + 1)
-    void adjacentNeighborsOf(const MultiIndex& position,
-                             std::array<std::pair<int, int>, space_dimension>& neighbors) const;
-    // Get diagonal neighbors of a node.
-    // If the neighbor index is -1, then there is no neighboring node in that direction.
-    // - The first entry (neighbors[0]) represents:
-    //   - First: bottom left neighbor (r - 1, theta - 1)
-    //   - Second: bottom right neighbor (r + 1, theta - 1)
-    // - The second entry (neighbors[1]) represents:
-    //   - First: top left neighbor (r - 1, theta + 1)
-    //   - Second: top right neighbor (r + 1, theta + 1)
-    void diagonalNeighborsOf(const MultiIndex& position,
-                             std::array<std::pair<int, int>, space_dimension>& neighbors) const;
-    // Neighbor distances
-    // Get distances to adjacent neighbors of a node.
-    // If there is no neighboring node in that direction, then the neighbor distance is 0.
-    void adjacentNeighborDistances(const MultiIndex& position,
-                                   std::array<std::pair<double, double>, space_dimension>& neighbor_distances) const;
-
 private:
-    // --------------- //
-    // Private members //
-    // --------------- //
-
-    // We will use the convention:
+    // We use the convention:
     // radii = [R0, ..., R], angles = [0, ..., 2*pi]
     // Note that ntheta will be one less than the size of angles since 0 and 2pi are the same point.
     int nr_; // number of nodes in radial direction
@@ -155,15 +104,11 @@ class PolarGrid
      * - numberCircularSmootherNodes + numberRadialSmootherNodes = number_of_nodes()
      */
 
-    // ------------------------ //
-    // Private Helper Functions //
-    // ------------------------ //
-
     // Check parameter validity
     void checkParameters(const std::vector<double>& radii, const std::vector<double>& angles) const;
+
     // Initialize radial_spacings_, angular_spacings_
     void initializeDistances();
-
     // Initializes line splitting parameters for Circle/radial indexing.
     // splitting_radius: The radius value used for dividing the smoother into a circular and radial section.
     //      If std::nullopt, automatic line-splitting is enabled.
@@ -192,4 +137,4 @@ class PolarGrid
 // ---------------------------------------------------- //
 PolarGrid coarseningGrid(const PolarGrid& grid);
 
-#include "polargrid.inl" // Include the inline function definitions
\ No newline at end of file
+#include "polargrid.inl" // Include the inline function definitions
diff --git a/include/PolarGrid/polargrid.inl b/include/PolarGrid/polargrid.inl
index 512343dc..9d9df122 100755
--- a/include/PolarGrid/polargrid.inl
+++ b/include/PolarGrid/polargrid.inl
@@ -28,6 +28,11 @@ inline double PolarGrid::theta(const int theta_index) const
     return angles_[theta_index];
 }
 
+inline double PolarGrid::smootherSplittingRadius() const
+{
+    return smoother_splitting_radius_;
+}
+
 // Get the number of circles in the circular smoother.
 inline int PolarGrid::numberSmootherCircles() const
 {
@@ -60,14 +65,9 @@ inline double PolarGrid::angularSpacing(const int unwrapped_theta_index) const
 {
     // unwrapped_theta_index may be negative or larger than ntheta() to allow for periodicity.
     const int theta_index = wrapThetaIndex(unwrapped_theta_index);
-    assert(theta_index >= 0 && theta_index < ntheta());
     return angular_spacings_[theta_index];
 }
 
-// ------------------ //
-// Optimized indexing //
-// ------------------ //
-
 inline int PolarGrid::wrapThetaIndex(const int unwrapped_theta_index) const
 {
     // The unwrapped_theta_index may be negative or exceed the number of theta steps (ntheta()),
@@ -80,40 +80,34 @@ inline int PolarGrid::wrapThetaIndex(const int unwrapped_theta_index) const
     //   This effectively computes unwrapped_theta_index % ntheta(), because it discards all higher bits.
     //
     // If ntheta is not a power of two, we use the standard modulo approach to handle wrapping.
-    return is_ntheta_PowerOfTwo_ ? unwrapped_theta_index & (ntheta() - 1) : (unwrapped_theta_index % ntheta() + ntheta()) % ntheta();
+    int theta_index = is_ntheta_PowerOfTwo_ ? unwrapped_theta_index & (ntheta() - 1)
+                                            : (unwrapped_theta_index % ntheta() + ntheta()) % ntheta();
+    assert(0 <= theta_index && theta_index < ntheta());
+    return theta_index;
 }
 
 inline int PolarGrid::index(const int r_index, const int unwrapped_theta_index) const
 {
     // unwrapped_theta_index may be negative or larger than ntheta() to allow for periodicity.
     assert(0 <= r_index && r_index < nr());
-    const int theta_index = wrapThetaIndex(unwrapped_theta_index);
-    assert(0 <= theta_index && theta_index < ntheta());
-    return r_index < numberSmootherCircles()
-               ? theta_index + ntheta() * r_index
-               : numberCircularSmootherNodes() + r_index - numberSmootherCircles() + lengthSmootherRadial() * theta_index;
-}
-
-inline int PolarGrid::fastIndex(const int r_index, const int theta_index) const
-{
-    assert(0 <= r_index && r_index < nr());
-    assert(0 <= theta_index && theta_index < ntheta());
-    return r_index < numberSmootherCircles()
-               ? theta_index + ntheta() * r_index
-               : numberCircularSmootherNodes() + r_index - numberSmootherCircles() + lengthSmootherRadial() * theta_index;
+    int theta_index = wrapThetaIndex(unwrapped_theta_index);
+    int global_index =
+        r_index < numberSmootherCircles()
+            ? theta_index + ntheta() * r_index
+            : numberCircularSmootherNodes() + r_index - numberSmootherCircles() + lengthSmootherRadial() * theta_index;
+    assert(0 <= global_index && global_index < numberOfNodes());
+    return global_index;
 }
 
 inline void PolarGrid::multiIndex(const int node_index, int& r_index, int& theta_index) const
 {
     assert(0 <= node_index && node_index < numberOfNodes());
-    if (node_index < numberCircularSmootherNodes())
-    {
-        r_index = node_index / ntheta();
-        theta_index = is_ntheta_PowerOfTwo_ ? node_index & (ntheta() - 1) : node_index % ntheta();
+    if (node_index < numberCircularSmootherNodes()) {
+        r_index     = node_index / ntheta();
+        theta_index = wrapThetaIndex(node_index);
     }
-    else
-    {
+    else {
         theta_index = (node_index - numberCircularSmootherNodes()) / lengthSmootherRadial();
-        r_index = numberSmootherCircles() + (node_index - numberCircularSmootherNodes()) % lengthSmootherRadial();
+        r_index     = numberSmootherCircles() + (node_index - numberCircularSmootherNodes()) % lengthSmootherRadial();
     }
 }
diff --git a/include/common/space_dimension.h b/include/common/space_dimension.h
deleted file mode 100644
index 9ba369e1..00000000
--- a/include/common/space_dimension.h
+++ /dev/null
@@ -1,3 +0,0 @@
-#pragma once
-
-constexpr inline int space_dimension = 2;
\ No newline at end of file
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 1b7ea608..eb06eece 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -6,10 +6,8 @@ add_subdirectory(InputFunctions)
 # Gather all source files
 # file(GLOB_RECURSE POLAR_GRID_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/*.cpp)
 set(POLAR_GRID_SOURCES
-    ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/anisotropic_division.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/multiindex.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/point.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/polargrid.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/PolarGrid/anisotropic_division.cpp
 )
 
 # file(GLOB_RECURSE GMG_POLAR_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/GMGPolar/*.cpp)
diff --git a/src/GMGPolar/setup.cpp b/src/GMGPolar/setup.cpp
index f331c14c..3c897a1b 100644
--- a/src/GMGPolar/setup.cpp
+++ b/src/GMGPolar/setup.cpp
@@ -271,7 +271,8 @@ void GMGPolar::printSettings() const
     const PolarGrid& finest_grid   = levels_.front().grid();
     const PolarGrid& coarsest_grid = levels_.back().grid();
 
-    std::cout << "r ∈ [" << finest_grid.radii().front() << ", " << finest_grid.radii().back() << "], θ ∈ [0, 2π]\n";
+    std::cout << "r ∈ [" << finest_grid.radius(0) << ", " << finest_grid.radius(finest_grid.nr() - 1)
+              << "], θ ∈ [0, 2π]\n";
     std::cout << "(nr × nθ) = (" << finest_grid.nr() << " × " << finest_grid.ntheta() << ") → (" << coarsest_grid.nr()
               << " × " << coarsest_grid.ntheta() << ")\n";
     std::cout << "Smoother: " << finest_grid.numberSmootherCircles() << " circles\n";
diff --git a/src/PolarGrid/multiindex.cpp b/src/PolarGrid/multiindex.cpp
deleted file mode 100644
index ec3c16a3..00000000
--- a/src/PolarGrid/multiindex.cpp
+++ /dev/null
@@ -1,71 +0,0 @@
-#include "../../include/PolarGrid/multiindex.h"
-
-MultiIndex::MultiIndex(int value)
-{
-    assert(space_dimension >= 0);
-    std::fill(std::begin(data_), std::end(data_), value);
-}
-
-MultiIndex::MultiIndex(int i, int j)
-{
-    assert(space_dimension == 2);
-    data_[0] = i;
-    data_[1] = j;
-}
-
-MultiIndex::MultiIndex(int i, int j, int k)
-{
-    assert(space_dimension == 3);
-    data_[0] = i;
-    data_[1] = j;
-    data_[2] = k;
-}
-
-int MultiIndex::size() const
-{
-    return space_dimension;
-}
-
-bool MultiIndex::operator==(const MultiIndex& other) const
-{
-    for (int d = 0; d < space_dimension; d++) {
-        if (data_[d] != other.data_[d]) {
-            return false;
-        }
-    }
-    return true;
-}
-
-bool MultiIndex::operator!=(const MultiIndex& other) const
-{
-    for (int d = 0; d < space_dimension; d++) {
-        if (data_[d] != other.data_[d]) {
-            return true;
-        }
-    }
-    return false;
-}
-
-const int& MultiIndex::operator[](int i) const
-{
-    assert(i >= 0);
-    assert(i < space_dimension);
-    return data_[i];
-}
-
-int& MultiIndex::operator[](int i)
-{
-    assert(i >= 0);
-    assert(i < space_dimension);
-    return data_[i];
-}
-
-int* MultiIndex::data()
-{
-    return data_;
-}
-
-const int* MultiIndex::data() const
-{
-    return data_;
-}
\ No newline at end of file
diff --git a/src/PolarGrid/point.cpp b/src/PolarGrid/point.cpp
deleted file mode 100644
index cb0bdf75..00000000
--- a/src/PolarGrid/point.cpp
+++ /dev/null
@@ -1,85 +0,0 @@
-#include "../../include/PolarGrid/point.h"
-
-Point::Point(double value)
-{
-    assert(space_dimension >= 0);
-    std::fill(std::begin(data_), std::end(data_), value);
-}
-
-Point::Point(double x, double y)
-{
-    assert(space_dimension == 2);
-    data_[0] = x;
-    data_[1] = y;
-}
-
-Point::Point(double x, double y, double z)
-{
-    assert(space_dimension == 3);
-    data_[0] = x;
-    data_[1] = y;
-    data_[2] = z;
-}
-
-int Point::size() const
-{
-    return space_dimension;
-}
-
-double Point::operator[](int i) const
-{
-    assert(i >= 0);
-    assert(i < space_dimension);
-    return data_[i];
-}
-
-double& Point::operator[](int i)
-{
-    assert(i >= 0);
-    assert(i < space_dimension);
-    return data_[i];
-}
-
-bool equals(const Point& lhs, const Point& rhs)
-{
-    for (int i = 0; i < space_dimension; ++i) {
-        if (!equals(lhs[i], rhs[i]))
-            return false;
-    }
-    return true;
-}
-
-double norm(const Point& point)
-{
-    double result = 0;
-    for (int i = 0; i < space_dimension; ++i) {
-        result += point[i] * point[i];
-    }
-    return sqrt(result);
-}
-
-void add(Point& result, const Point& lhs, const Point& rhs)
-{
-    assert(lhs.size() == rhs.size());
-    assert(result.size() == lhs.size());
-    for (int i = 0; i < space_dimension; i++) {
-        result[i] = lhs[i] + rhs[i];
-    }
-}
-
-void subtract(Point& result, const Point& lhs, const Point& rhs)
-{
-    assert(lhs.size() == rhs.size());
-    assert(result.size() == lhs.size());
-    for (int i = 0; i < space_dimension; i++) {
-        result[i] = lhs[i] - rhs[i];
-    }
-}
-
-void multiply(Point& result, double scalar, const Point& lhs)
-{
-    assert(result.size() == lhs.size());
-    for (int i = 0; i < space_dimension; i++) {
-        result[i] = scalar * lhs[i];
-    }
-}
\ No newline at end of file
diff --git a/src/PolarGrid/polargrid.cpp b/src/PolarGrid/polargrid.cpp
index 2059640d..5f658aec 100644
--- a/src/PolarGrid/polargrid.cpp
+++ b/src/PolarGrid/polargrid.cpp
@@ -235,143 +235,9 @@ PolarGrid coarseningGrid(const PolarGrid& fineGrid)
     }
 }
 
-// ---------------- //
-// Getter Functions //
-// ---------------- //
-
-const std::vector<double>& PolarGrid::radii() const
-{
-    return radii_;
-}
-const std::vector<double>& PolarGrid::angles() const
-{
-    return angles_;
-}
-
-// Get the radius at which the grid is split into circular and radial smoothing
-double PolarGrid::smootherSplittingRadius() const
-{
-    return smoother_splitting_radius_;
-}
-
-// ------------------------------------- //
-// Definition of node indexing.          //
-// Based on the circular-radial smoother //
-// ------------------------------------- //
-
-/* OPTIMIZED INDEXING IS DEFINED IN include/PolarGrid/polargrid.inl */
-
-int PolarGrid::index(const MultiIndex& position) const
-{
-    assert(position[0] >= 0 && position[0] < nr());
-    assert(position[1] >= 0 && position[1] < ntheta());
-    if (position[0] < numberSmootherCircles()) {
-        return position[1] + ntheta() * position[0];
-    }
-    else {
-        return numberCircularSmootherNodes() +
-               (position[0] - numberSmootherCircles() + lengthSmootherRadial() * position[1]);
-    }
-}
-
-MultiIndex PolarGrid::multiIndex(const int node_index) const
-{
-    assert(0 <= node_index && node_index < numberOfNodes());
-    if (node_index < numberCircularSmootherNodes()) {
-        auto result = std::div(node_index, ntheta());
-        return MultiIndex(result.quot, result.rem);
-    }
-    else {
-        auto result = std::div(node_index - numberCircularSmootherNodes(), lengthSmootherRadial());
-        return MultiIndex(numberSmootherCircles() + result.rem, result.quot);
-    }
-}
-
-Point PolarGrid::polarCoordinates(const MultiIndex& position) const
-{
-    assert(position[0] >= 0 && position[0] < nr());
-    assert(position[1] >= 0 && position[1] < ntheta());
-    return Point(radii_[position[0]], angles_[(position[1])]);
-}
-
-void PolarGrid::adjacentNeighborDistances(
-    const MultiIndex& position, std::array<std::pair<double, double>, space_dimension>& neighbor_distance) const
-{
-    assert(position[0] >= 0 && position[0] < nr());
-    assert(position[1] >= 0 && position[1] < ntheta());
-
-    neighbor_distance[0].first  = (position[0] <= 0) ? 0.0 : radialSpacing(position[0] - 1);
-    neighbor_distance[0].second = (position[0] >= nr() - 1) ? 0.0 : radialSpacing(position[0]);
-
-    neighbor_distance[1].first  = angularSpacing(position[1] - 1);
-    neighbor_distance[1].second = angularSpacing(position[1]);
-}
-
-void PolarGrid::adjacentNeighborsOf(const MultiIndex& position,
-                                    std::array<std::pair<int, int>, space_dimension>& neighbors) const
-{
-    assert(position[0] >= 0 && position[0] < nr());
-    assert(position[1] >= 0 && position[1] < ntheta());
-
-    MultiIndex neigbor_position = position;
-    neigbor_position[0] -= 1;
-    neighbors[0].first = (neigbor_position[0] < 0) ? -1 : index(neigbor_position);
-
-    neigbor_position = position;
-    neigbor_position[0] += 1;
-    neighbors[0].second = (neigbor_position[0] >= nr()) ? -1 : index(neigbor_position);
-
-    neigbor_position = position;
-    neigbor_position[1] -= 1;
-    if (neigbor_position[1] < 0)
-        neigbor_position[1] += ntheta();
-    neighbors[1].first = index(neigbor_position);
-
-    neigbor_position = position;
-    neigbor_position[1] += 1;
-    if (neigbor_position[1] >= ntheta())
-        neigbor_position[1] -= ntheta();
-    neighbors[1].second = index(neigbor_position);
-}
-
-void PolarGrid::diagonalNeighborsOf(const MultiIndex& position,
-                                    std::array<std::pair<int, int>, space_dimension>& neighbors) const
-{
-    assert(position[0] >= 0 && position[0] < nr());
-    assert(position[1] >= 0 && position[1] < ntheta());
-
-    MultiIndex neigbor_position = position;
-    neigbor_position[0] -= 1;
-    neigbor_position[1] -= 1;
-    if (neigbor_position[1] < 0)
-        neigbor_position[1] += ntheta();
-    neighbors[0].first = (neigbor_position[0] < 0) ? -1 : index(neigbor_position);
-
-    neigbor_position = position;
-    neigbor_position[0] += 1;
-    neigbor_position[1] -= 1;
-    if (neigbor_position[1] < 0)
-        neigbor_position[1] += ntheta();
-    neighbors[0].second = (neigbor_position[0] >= nr()) ? -1 : index(neigbor_position);
-
-    neigbor_position = position;
-    neigbor_position[0] -= 1;
-    neigbor_position[1] += 1;
-    if (neigbor_position[1] >= ntheta())
-        neigbor_position[1] -= ntheta();
-    neighbors[1].first = (neigbor_position[0] < 0) ? -1 : index(neigbor_position);
-
-    neigbor_position = position;
-    neigbor_position[0] += 1;
-    neigbor_position[1] += 1;
-    if (neigbor_position[1] >= ntheta())
-        neigbor_position[1] -= ntheta();
-    neighbors[1].second = (neigbor_position[0] >= nr()) ? -1 : index(neigbor_position);
-}
-
 // ------------------------ //
 // Check parameter validity //
-// ---------------------..- //
+// ------------------------ //
 
 void PolarGrid::checkParameters(const std::vector<double>& radii, const std::vector<double>& angles) const
 {
diff --git a/tests/DirectSolver/directSolver.cpp b/tests/DirectSolver/directSolver.cpp
index 06e7e5ff..3505af85 100644
--- a/tests/DirectSolver/directSolver.cpp
+++ b/tests/DirectSolver/directSolver.cpp
@@ -98,8 +98,9 @@ TEST(DirectSolverTest, directSolver_DirBC_Interior)
 
     ASSERT_EQ(solution_Give.size(), solution_Take.size());
     for (uint index = 0; index < solution_Give.size(); index++) {
-        MultiIndex alpha = level.grid().multiIndex(index);
-        if (alpha[0] == 0 && !DirBC_Interior)
+        int i_r, i_theta;
+        level.grid().multiIndex(index, i_r, i_theta);
+        if (i_r == 0 && !DirBC_Interior)
             ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-11);
         else
             ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-11);
@@ -153,8 +154,9 @@ TEST(DirectSolverTest, directSolver_AcrossOrigin)
 
     ASSERT_EQ(solution_Give.size(), solution_Take.size());
     for (uint index = 0; index < solution_Give.size(); index++) {
-        MultiIndex alpha = level.grid().multiIndex(index);
-        if (alpha[0] == 0 && !DirBC_Interior)
+        int i_r, i_theta;
+        level.grid().multiIndex(index, i_r, i_theta);
+        if (i_r == 0 && !DirBC_Interior)
             ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-8);
         else
             ASSERT_NEAR(solution_Give(index), solution_Take(index), 1e-8);
diff --git a/tests/DirectSolver/directSolverNoMumps.cpp b/tests/DirectSolver/directSolverNoMumps.cpp
index 7ae3eec7..49274ad1 100644
--- a/tests/DirectSolver/directSolverNoMumps.cpp
+++ b/tests/DirectSolver/directSolverNoMumps.cpp
@@ -96,8 +96,9 @@ TEST(DirectSolverTestNoMumps, directSolver_DirBC_Interior)
 
     ASSERT_EQ(solution_Give.size(), solution_Take.size());
     for (uint index = 0; index < solution_Give.size(); index++) {
-        MultiIndex alpha = level.grid().multiIndex(index);
-        if (alpha[0] == 0 && !DirBC_Interior)
+        int i_r, i_theta;
+        level.grid().multiIndex(index, i_r, i_theta);
+        if (i_r == 0 && !DirBC_Interior)
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11);
         else
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11);
@@ -151,8 +152,9 @@ TEST(DirectSolverTestNoMumps, directSolver_AcrossOrigin)
 
     ASSERT_EQ(solution_Give.size(), solution_Take.size());
     for (uint index = 0; index < solution_Give.size(); index++) {
-        MultiIndex alpha = level.grid().multiIndex(index);
-        if (alpha[0] == 0 && !DirBC_Interior)
+        int i_r, i_theta;
+        level.grid().multiIndex(index, i_r, i_theta);
+        if (i_r == 0 && !DirBC_Interior)
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-8);
         else
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-9);
diff --git a/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp b/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp
index 714403a5..8af5fc06 100644
--- a/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp
+++ b/tests/ExtrapolatedSmoother/extrapolated_smoother.cpp
@@ -83,8 +83,9 @@ TEST(ExtrapolatedSmootherTest, extrapolatedSmoother_DirBC_Interior)
 
     ASSERT_EQ(solution_Give.size(), solution_Take.size());
     for (uint index = 0; index < solution_Give.size(); index++) {
-        MultiIndex alpha = level.grid().multiIndex(index);
-        if (alpha[0] == 0 && !DirBC_Interior)
+        int i_r, i_theta;
+        level.grid().multiIndex(index, i_r, i_theta);
+        if (i_r == 0 && !DirBC_Interior)
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11);
         else
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11);
@@ -142,8 +143,9 @@ TEST(ExtrapolatedSmootherTest, extrapolatedSmoother_AcossOrigin)
 
     ASSERT_EQ(solution_Give.size(), solution_Take.size());
     for (uint index = 0; index < solution_Give.size(); index++) {
-        MultiIndex alpha = level.grid().multiIndex(index);
-        if (alpha[0] == 0 && !DirBC_Interior)
+        int i_r, i_theta;
+        level.grid().multiIndex(index, i_r, i_theta);
+        if (i_r == 0 && !DirBC_Interior)
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-8);
         else
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11);
diff --git a/tests/PolarGrid/polargrid.cpp b/tests/PolarGrid/polargrid.cpp
index c3915acc..31a87f2a 100644
--- a/tests/PolarGrid/polargrid.cpp
+++ b/tests/PolarGrid/polargrid.cpp
@@ -58,17 +58,18 @@ TEST(PolarGridTest, IndexingTest)
 
     for (int i = 0; i < grid.nr(); i++) {
         for (int j = 0; j < grid.ntheta(); j++) {
-            MultiIndex alpha(i, j);
-            int node_index  = grid.index(alpha);
-            MultiIndex beta = grid.multiIndex(node_index);
-            ASSERT_EQ(alpha[0], beta[0]);
-            ASSERT_EQ(alpha[1], beta[1]);
+            int node_index = grid.index(i, j);
+            int r_out, theta_out;
+            grid.multiIndex(node_index, r_out, theta_out);
+            ASSERT_EQ(i, r_out);
+            ASSERT_EQ(j, theta_out);
         }
     }
 
     for (int i = 0; i < grid.numberOfNodes(); i++) {
-        MultiIndex alpha = grid.multiIndex(i);
-        int node_index   = grid.index(alpha);
+        int r_out, theta_out;
+        grid.multiIndex(i, r_out, theta_out);
+        int node_index = grid.index(r_out, theta_out);
         ASSERT_EQ(node_index, i);
     }
 }
@@ -82,39 +83,39 @@ TEST(PolarGridTest, IndexingValuesTest)
     PolarGrid grid(radii, angles, splitting_radius);
 
     {
-        MultiIndex alpha(2, 6);
-        int node_index = grid.index(alpha);
+        int node_index = grid.index(2, 6);
         ASSERT_EQ(node_index, 22);
-        MultiIndex beta = grid.multiIndex(node_index);
-        ASSERT_EQ(alpha[0], beta[0]);
-        ASSERT_EQ(alpha[1], beta[1]);
+        int r_out, theta_out;
+        grid.multiIndex(node_index, r_out, theta_out);
+        ASSERT_EQ(2, r_out);
+        ASSERT_EQ(6, theta_out);
     }
 
     {
-        MultiIndex alpha(3, 2);
-        int node_index = grid.index(alpha);
+        int node_index = grid.index(3, 2);
         ASSERT_EQ(node_index, 26);
-        MultiIndex beta = grid.multiIndex(node_index);
-        ASSERT_EQ(alpha[0], beta[0]);
-        ASSERT_EQ(alpha[1], beta[1]);
+        int r_out, theta_out;
+        grid.multiIndex(node_index, r_out, theta_out);
+        ASSERT_EQ(3, r_out);
+        ASSERT_EQ(2, theta_out);
     }
 
     {
-        MultiIndex alpha(6, 4);
-        int node_index = grid.index(alpha);
+        int node_index = grid.index(6, 4);
         ASSERT_EQ(node_index, 54);
-        MultiIndex beta = grid.multiIndex(node_index);
-        ASSERT_EQ(alpha[0], beta[0]);
-        ASSERT_EQ(alpha[1], beta[1]);
+        int r_out, theta_out;
+        grid.multiIndex(node_index, r_out, theta_out);
+        ASSERT_EQ(6, r_out);
+        ASSERT_EQ(4, theta_out);
     }
 
     {
-        MultiIndex alpha(4, 7);
-        int node_index = grid.index(alpha);
+        int node_index = grid.index(4, 7);
         ASSERT_EQ(node_index, 67);
-        MultiIndex beta = grid.multiIndex(node_index);
-        ASSERT_EQ(alpha[0], beta[0]);
-        ASSERT_EQ(alpha[1], beta[1]);
+        int r_out, theta_out;
+        grid.multiIndex(node_index, r_out, theta_out);
+        ASSERT_EQ(4, r_out);
+        ASSERT_EQ(7, theta_out);
     }
 }
 
@@ -126,129 +127,27 @@ TEST(PolarGridTest, CoordinatesTest)
     double splitting_radius = 0.6;
     PolarGrid grid(radii, angles, splitting_radius);
 
-    MultiIndex alpha(3, 2);
-    Point P1 = grid.polarCoordinates(alpha);
-    ASSERT_DOUBLE_EQ(P1[0], 0.5);
-    ASSERT_DOUBLE_EQ(P1[1], M_PI / 8);
+    ASSERT_DOUBLE_EQ(grid.radius(3), 0.5);
+    ASSERT_DOUBLE_EQ(grid.theta(2), M_PI / 8);
 
-    alpha[0] += 4;
-    alpha[1] -= 1;
-    P1 = grid.polarCoordinates(alpha);
-    ASSERT_DOUBLE_EQ(P1[0], 1.4);
-    ASSERT_DOUBLE_EQ(P1[1], M_PI / 16);
+    ASSERT_DOUBLE_EQ(grid.radius(7), 1.4);
+    ASSERT_DOUBLE_EQ(grid.theta(1), M_PI / 16);
 }
 
-TEST(PolarGridTest, NeighborDistanceTest)
+TEST(PolarGridTest, SpacingTest)
 {
     std::vector<double> radii  = {0.1, 0.2, 0.25, 0.5, 0.8, 0.9, 1.3, 1.4, 2.0};
     std::vector<double> angles = {
         0, M_PI / 16, M_PI / 8, M_PI / 2, M_PI, M_PI + M_PI / 16, M_PI + M_PI / 8, M_PI + M_PI / 2, M_PI + M_PI};
-    double splitting_radius = 0.6;
-    PolarGrid grid(radii, angles, splitting_radius);
-
-    std::array<std::pair<double, double>, space_dimension> neighbor_distance;
-    {
-        MultiIndex alpha(3, 2);
-        grid.adjacentNeighborDistances(alpha, neighbor_distance);
-        ASSERT_DOUBLE_EQ(neighbor_distance[0].first, 0.5 - 0.25);
-        ASSERT_DOUBLE_EQ(neighbor_distance[0].second, 0.8 - 0.5);
-        ASSERT_DOUBLE_EQ(neighbor_distance[1].first, M_PI / 8 - M_PI / 16);
-        ASSERT_DOUBLE_EQ(neighbor_distance[1].second, M_PI / 2 - M_PI / 8);
-    }
-
-    {
-        MultiIndex alpha(8, 7);
-        grid.adjacentNeighborDistances(alpha, neighbor_distance);
-        ASSERT_DOUBLE_EQ(neighbor_distance[0].first, 2.0 - 1.4);
-        ASSERT_DOUBLE_EQ(neighbor_distance[0].second, 0.0);
-        ASSERT_DOUBLE_EQ(neighbor_distance[1].first, M_PI / 2 - M_PI / 8);
-        ASSERT_DOUBLE_EQ(neighbor_distance[1].second, M_PI / 2);
-    }
-
-    {
-        MultiIndex alpha(4, 0);
-        grid.adjacentNeighborDistances(alpha, neighbor_distance);
-        ASSERT_DOUBLE_EQ(neighbor_distance[0].first, 0.8 - 0.5);
-        ASSERT_DOUBLE_EQ(neighbor_distance[0].second, 0.9 - 0.8);
-        ASSERT_DOUBLE_EQ(neighbor_distance[1].first, M_PI / 2);
-        ASSERT_DOUBLE_EQ(neighbor_distance[1].second, M_PI / 16);
-    }
-
-    {
-        MultiIndex alpha(0, 5);
-        grid.adjacentNeighborDistances(alpha, neighbor_distance);
-        ASSERT_DOUBLE_EQ(neighbor_distance[0].first, 0.0);
-        ASSERT_DOUBLE_EQ(neighbor_distance[0].second, 0.2 - 0.1);
-        ASSERT_NEAR(neighbor_distance[1].first, M_PI / 16, 1e-15);
-        ASSERT_NEAR(neighbor_distance[1].second, M_PI / 8 - M_PI / 16, 1e-15);
-        ASSERT_EQ(equals(neighbor_distance[1].first, M_PI / 16), true);
-        ASSERT_EQ(equals(neighbor_distance[1].second, M_PI / 8 - M_PI / 16), true);
-    }
-}
-
-TEST(PolarGridTest, NeighborsTest)
-{
-    std::vector<double> radii  = {0.1, 0.2, 0.25, 0.5, 0.8, 0.9, 1.3, 1.4, 2.0};
-    std::vector<double> angles = {
-        0, M_PI / 16, M_PI / 8, M_PI / 2, M_PI, M_PI + M_PI / 16, M_PI + M_PI / 8, M_PI + M_PI / 2, M_PI + M_PI};
-    double splitting_radius = 0.6;
-    PolarGrid grid(radii, angles, splitting_radius);
-
-    std::array<std::pair<int, int>, space_dimension> neighbors;
-
-    {
-        MultiIndex alpha(3, 2);
-        grid.adjacentNeighborsOf(alpha, neighbors);
-        ASSERT_EQ(neighbors[0].first, 18);
-        ASSERT_EQ(neighbors[0].second, 42);
-        ASSERT_EQ(neighbors[1].first, 25);
-        ASSERT_EQ(neighbors[1].second, 27);
-        grid.diagonalNeighborsOf(alpha, neighbors);
-        ASSERT_EQ(neighbors[0].first, 17);
-        ASSERT_EQ(neighbors[0].second, 37);
-        ASSERT_EQ(neighbors[1].first, 19);
-        ASSERT_EQ(neighbors[1].second, 47);
-    }
-
-    {
-        MultiIndex alpha(8, 7);
-        grid.adjacentNeighborsOf(alpha, neighbors);
-        ASSERT_EQ(neighbors[0].first, 70);
-        ASSERT_EQ(neighbors[0].second, -1);
-        ASSERT_EQ(neighbors[1].first, 66);
-        ASSERT_EQ(neighbors[1].second, 36);
-        grid.diagonalNeighborsOf(alpha, neighbors);
-        ASSERT_EQ(neighbors[0].first, 65);
-        ASSERT_EQ(neighbors[0].second, -1);
-        ASSERT_EQ(neighbors[1].first, 35);
-        ASSERT_EQ(neighbors[1].second, -1);
-    }
+    PolarGrid grid(radii, angles);
 
-    {
-        MultiIndex alpha(4, 0);
-        grid.adjacentNeighborsOf(alpha, neighbors);
-        ASSERT_EQ(neighbors[0].first, 24);
-        ASSERT_EQ(neighbors[0].second, 33);
-        ASSERT_EQ(neighbors[1].first, 67);
-        ASSERT_EQ(neighbors[1].second, 37);
-        grid.diagonalNeighborsOf(alpha, neighbors);
-        ASSERT_EQ(neighbors[0].first, 31);
-        ASSERT_EQ(neighbors[0].second, 68);
-        ASSERT_EQ(neighbors[1].first, 25);
-        ASSERT_EQ(neighbors[1].second, 38);
-    }
+    // Test radial spacings
+    ASSERT_DOUBLE_EQ(grid.radialSpacing(0), 0.2 - 0.1);
+    ASSERT_DOUBLE_EQ(grid.radialSpacing(1), 0.25 - 0.2);
+    ASSERT_DOUBLE_EQ(grid.radialSpacing(2), 0.5 - 0.25);
 
-    {
-        MultiIndex alpha(0, 5);
-        grid.adjacentNeighborsOf(alpha, neighbors);
-        ASSERT_EQ(neighbors[0].first, -1);
-        ASSERT_EQ(neighbors[0].second, 13);
-        ASSERT_EQ(neighbors[1].first, 4);
-        ASSERT_EQ(neighbors[1].second, 6);
-        grid.diagonalNeighborsOf(alpha, neighbors);
-        ASSERT_EQ(neighbors[0].first, -1);
-        ASSERT_EQ(neighbors[0].second, 12);
-        ASSERT_EQ(neighbors[1].first, -1);
-        ASSERT_EQ(neighbors[1].second, 14);
-    }
+    // Test angular spacings
+    ASSERT_DOUBLE_EQ(grid.angularSpacing(0), M_PI / 16 - 0);
+    ASSERT_DOUBLE_EQ(grid.angularSpacing(1), M_PI / 8 - M_PI / 16);
+    ASSERT_DOUBLE_EQ(grid.angularSpacing(2), M_PI / 2 - M_PI / 8);
 }
diff --git a/tests/Residual/residual.cpp b/tests/Residual/residual.cpp
index 70ec6449..4bc72bb1 100644
--- a/tests/Residual/residual.cpp
+++ b/tests/Residual/residual.cpp
@@ -73,8 +73,9 @@ TEST(OperatorATest, applyA_DirBC_Interior)
 
     ASSERT_EQ(result_Give.size(), result_Take.size());
     for (uint index = 0; index < result_Give.size(); index++) {
-        MultiIndex alpha = level.grid().multiIndex(index);
-        if (alpha[0] == 0 && !DirBC_Interior)
+        int i_r, i_theta;
+        level.grid().multiIndex(index, i_r, i_theta);
+        if (i_r == 0 && !DirBC_Interior)
             ASSERT_NEAR(result_Give[index], result_Take[index], 1e-8);
         else
             ASSERT_NEAR(result_Give[index], result_Take[index], 1e-11);
@@ -129,8 +130,9 @@ TEST(OperatorATest, applyA_AcrossOrigin)
 
     ASSERT_EQ(result_Give.size(), result_Take.size());
     for (uint index = 0; index < result_Give.size(); index++) {
-        MultiIndex alpha = level.grid().multiIndex(index);
-        if (alpha[0] == 0 && !DirBC_Interior)
+        int i_r, i_theta;
+        level.grid().multiIndex(index, i_r, i_theta);
+        if (i_r == 0 && !DirBC_Interior)
             ASSERT_NEAR(result_Give[index], result_Take[index], 1e-8);
         else
             ASSERT_NEAR(result_Give[index], result_Take[index], 1e-11);
diff --git a/tests/Smoother/smoother.cpp b/tests/Smoother/smoother.cpp
index 8fad3f54..9e6cc214 100644
--- a/tests/Smoother/smoother.cpp
+++ b/tests/Smoother/smoother.cpp
@@ -80,8 +80,9 @@ TEST(SmootherTest, smoother_DirBC_Interior)
 
     ASSERT_EQ(solution_Give.size(), solution_Take.size());
     for (uint index = 0; index < solution_Give.size(); index++) {
-        MultiIndex alpha = level.grid().multiIndex(index);
-        if (alpha[0] == 0 && !DirBC_Interior)
+        int i_r, i_theta;
+        level.grid().multiIndex(index, i_r, i_theta);
+        if (i_r == 0 && !DirBC_Interior)
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11);
         else
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-11);
@@ -139,8 +140,9 @@ TEST(SmootherTest, smoother_AcrossOrigin)
 
     ASSERT_EQ(solution_Give.size(), solution_Take.size());
     for (uint index = 0; index < solution_Give.size(); index++) {
-        MultiIndex alpha = level.grid().multiIndex(index);
-        if (alpha[0] == 0 && !DirBC_Interior)
+        int i_r, i_theta;
+        level.grid().multiIndex(index, i_r, i_theta);
+        if (i_r == 0 && !DirBC_Interior)
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-8);
         else
             ASSERT_NEAR(solution_Give[index], solution_Take[index], 1e-10);