diff --git a/Applications/DataExplorer/mainwindow.cpp b/Applications/DataExplorer/mainwindow.cpp
index 4568081509aeddea4a46c2db9757c5d1e853bcac..591fbc0657f4843d6a8515375c87f6b3e973fd3e 100644
--- a/Applications/DataExplorer/mainwindow.cpp
+++ b/Applications/DataExplorer/mainwindow.cpp
@@ -1335,7 +1335,8 @@ void MainWindow::showMeshAnalysisDialog()
 void MainWindow::convertPointsToStations(std::string const& geo_name)
 {
     std::string stn_name = geo_name + " Stations";
-    int ret = _project.getGEOObjects().geoPointsToStations(geo_name, stn_name);
+    int ret = GeoLib::geoPointsToStations(_project.getGEOObjects(), geo_name,
+                                          stn_name);
     if (ret == 1)
     {
         OGSError::box("No points found to convert.");
diff --git a/Applications/FileIO/Gmsh/GMSHInterface.cpp b/Applications/FileIO/Gmsh/GMSHInterface.cpp
index 6b70c66057f91054a7f4246c091b36c6ef8acbcd..d49977821c12a4cada6390c99e7f50b36805d127 100644
--- a/Applications/FileIO/Gmsh/GMSHInterface.cpp
+++ b/Applications/FileIO/Gmsh/GMSHInterface.cpp
@@ -25,6 +25,7 @@
 #include "GeoLib/Polygon.h"
 #include "GeoLib/PolygonWithSegmentMarker.h"
 #include "GeoLib/PolylineWithSegmentMarker.h"
+#include "GeoLib/Station.h"
 #include "InfoLib/GitInfo.h"
 
 namespace FileIO
diff --git a/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp b/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
index 6d4a03c928bc83bd3ad58b95a39d6d0cc2419d19..b66547329fd8789f7575acf208e09f2613c65377 100644
--- a/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
+++ b/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
@@ -18,6 +18,7 @@
 #include "GeoLib/Polygon.h"
 #include "GeoLib/PolygonWithSegmentMarker.h"
 #include "GeoLib/PolylineWithSegmentMarker.h"
+#include "GeoLib/Station.h"
 
 namespace FileIO
 {
diff --git a/Applications/FileIO/Legacy/OGSIOVer4.cpp b/Applications/FileIO/Legacy/OGSIOVer4.cpp
index b8dc94ba004634d4a46ba19ab5094a42272030e2..9196c54843409487d09733d721161cb0ddae0cb2 100644
--- a/Applications/FileIO/Legacy/OGSIOVer4.cpp
+++ b/Applications/FileIO/Legacy/OGSIOVer4.cpp
@@ -29,6 +29,7 @@
 #include "GeoLib/PointVec.h"
 #include "GeoLib/Polygon.h"
 #include "GeoLib/Polyline.h"
+#include "GeoLib/Station.h"
 #include "GeoLib/Surface.h"
 #include "GeoLib/Triangle.h"
 
diff --git a/GeoLib/ClosestPair.h b/GeoLib/ClosestPair.h
deleted file mode 100644
index bf3835e91a6419c9fda3ef08eb776a09839ff25a..0000000000000000000000000000000000000000
--- a/GeoLib/ClosestPair.h
+++ /dev/null
@@ -1,38 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date   2011-01-25
- * \brief  Definition of the ClosestPair class.
- *
- * \copyright
- * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-// STL
-#include <vector>
-
-// GeoLib
-#include "Point.h"
-
-namespace GeoLib {
-
-class ClosestPair
-{
-public:
-    ClosestPair (std::vector<GeoLib::Point*> const & pnts, std::size_t id0, std::size_t id1) :
-        _pnts (pnts), _id0 (id0), _id1 (id1)
-    {}
-
-protected:
-    std::vector<GeoLib::Point*> const & _pnts;
-    std::size_t _id0;
-    std::size_t _id1;
-};
-
-} // end namespace GeoLib
diff --git a/GeoLib/GEOObjects.cpp b/GeoLib/GEOObjects.cpp
index d831f6b1a05ac6b839062fe439ceaaa97f7a8020..0465821bcc93c99a2d917b6e26f23ce08711dada 100644
--- a/GeoLib/GEOObjects.cpp
+++ b/GeoLib/GEOObjects.cpp
@@ -18,10 +18,15 @@
 
 #include "BaseLib/Logging.h"
 #include "BaseLib/StringTools.h"
+#include "Station.h"
 #include "Triangle.h"
 
 namespace GeoLib
 {
+void markUnusedPoints(GEOObjects const& geo_objects,
+                      std::string const& geo_name,
+                      std::vector<bool>& transfer_pnts);
+
 GEOObjects::GEOObjects() = default;
 
 GEOObjects::~GEOObjects()
@@ -650,10 +655,12 @@ void GEOObjects::renameGeometry(std::string const& old_name,
     }
 }
 
-void GEOObjects::markUnusedPoints(std::string const& geo_name,
-                                  std::vector<bool>& transfer_pnts) const
+void markUnusedPoints(GEOObjects const& geo_objects,
+                      std::string const& geo_name,
+                      std::vector<bool>& transfer_pnts)
 {
-    GeoLib::PolylineVec const* const ply_obj(getPolylineVecObj(geo_name));
+    GeoLib::PolylineVec const* const ply_obj(
+        geo_objects.getPolylineVecObj(geo_name));
     if (ply_obj)
     {
         std::vector<GeoLib::Polyline*> const& lines(*ply_obj->getVector());
@@ -667,7 +674,8 @@ void GEOObjects::markUnusedPoints(std::string const& geo_name,
         }
     }
 
-    GeoLib::SurfaceVec const* const sfc_obj(getSurfaceVecObj(geo_name));
+    GeoLib::SurfaceVec const* const sfc_obj(
+        geo_objects.getSurfaceVecObj(geo_name));
     if (sfc_obj)
     {
         std::vector<GeoLib::Surface*> const& surfaces = *sfc_obj->getVector();
@@ -685,11 +693,10 @@ void GEOObjects::markUnusedPoints(std::string const& geo_name,
     }
 }
 
-int GEOObjects::geoPointsToStations(std::string const& geo_name,
-                                    std::string& stn_name,
-                                    bool const only_unused_pnts)
+int geoPointsToStations(GEOObjects& geo_objects, std::string const& geo_name,
+                        std::string& stn_name, bool const only_unused_pnts)
 {
-    GeoLib::PointVec const* const pnt_obj(getPointVecObj(geo_name));
+    GeoLib::PointVec const* const pnt_obj(geo_objects.getPointVecObj(geo_name));
     if (pnt_obj == nullptr)
     {
         ERR("Point vector {:s} not found.", geo_name);
@@ -705,7 +712,7 @@ int GEOObjects::geoPointsToStations(std::string const& geo_name,
     std::vector<bool> transfer_pnts(n_pnts, true);
     if (only_unused_pnts)
     {
-        markUnusedPoints(geo_name, transfer_pnts);
+        markUnusedPoints(geo_objects, geo_name, transfer_pnts);
     }
 
     auto stations = std::make_unique<std::vector<GeoLib::Point*>>();
@@ -724,7 +731,7 @@ int GEOObjects::geoPointsToStations(std::string const& geo_name,
     }
     if (!stations->empty())
     {
-        addStationVec(std::move(stations), stn_name);
+        geo_objects.addStationVec(std::move(stations), stn_name);
     }
     else
     {
diff --git a/GeoLib/GEOObjects.h b/GeoLib/GEOObjects.h
index 4f90ec9ec07636adb1704fef8564778bf51e95b3..ae868285fc5e94bb6c5b21a6a3f069736d8afd53 100644
--- a/GeoLib/GEOObjects.h
+++ b/GeoLib/GEOObjects.h
@@ -25,9 +25,6 @@
 #include "PolylineVec.h"
 #include "Surface.h"
 #include "SurfaceVec.h"
-
-#include "Station.h"
-
 #include "GeoType.h"
 
 namespace GeoLib
@@ -256,15 +253,6 @@ public:
     void renameGeometry(std::string const& old_name,
                         std::string const& new_name);
 
-    /// Constructs a station-vector based on the points of a given geometry.
-    // @param geo_name name of the geometry
-    // @param stn_name name of the new station vector
-    // @param only_usused_pnts if true only points not in a line or surface are
-    // transferred, otherwise all points
-    int geoPointsToStations(std::string const& geo_name,
-                            std::string& stn_name,
-                            bool const only_unused_pnts = true);
-
     /// Returns the geo object for a geometric item of the given name and type for the associated geometry.
     const GeoLib::GeoObject* getGeoObject(const std::string &geo_name,
                                           GeoLib::GEOTYPE type,
@@ -374,8 +362,16 @@ private:
     void mergeSurfaces(std::vector<std::string> const& geo_names,
                        std::string const& merged_geo_name,
                        std::vector<std::size_t> const& pnt_offsets);
-
-    void markUnusedPoints(std::string const& geo_name,
-                          std::vector<bool>& transfer_pnts) const;
 };
+
+/// Constructs a station-vector based on the points of a given geometry.
+// @param geo_obj geometry manager object
+// @param geo_name name of the geometry
+// @param stn_name name of the new station vector
+// @param only_usused_pnts if true only points not in a line or surface are
+// transferred, otherwise all points
+int geoPointsToStations(GEOObjects& geo_objects, std::string const& geo_name,
+                        std::string& stn_name,
+                        bool const only_unused_pnts = true);
+
 }  // namespace GeoLib
diff --git a/GeoLib/Grid.h b/GeoLib/Grid.h
index 0d6900a6df8654248c8f0602a8ba5ba05ee08a70..442e63302f7000fd705604abd045ca68a0f163e4 100644
--- a/GeoLib/Grid.h
+++ b/GeoLib/Grid.h
@@ -95,10 +95,9 @@ public:
     getPntVecsOfGridCellsIntersectingCube(P const& center,
                                           double half_len) const;
 
-    void getPntVecsOfGridCellsIntersectingCuboid(
-        Eigen::Vector3d const& min_pnt,
-        Eigen::Vector3d const& max_pnt,
-        std::vector<std::vector<POINT*> const*>& pnts) const;
+    std::vector<std::vector<POINT*> const*>
+    getPntVecsOfGridCellsIntersectingCuboid(
+        Eigen::Vector3d const& min_pnt, Eigen::Vector3d const& max_pnt) const;
 
 #ifndef NDEBUG
     /**
@@ -253,30 +252,27 @@ std::vector<std::vector<POINT*> const*>
 Grid<POINT>::getPntVecsOfGridCellsIntersectingCube(P const& center,
                                                    double half_len) const
 {
-    std::vector<std::vector<POINT*> const*> pnts;
-    POINT tmp_pnt{center[0] - half_len, center[1] - half_len,
-                  center[2] - half_len};  // min
-    std::array<std::size_t, 3> min_coords(getGridCoords(tmp_pnt));
+    Eigen::Vector3d const c{center[0], center[1], center[2]};
+    std::array<std::size_t, 3> min_coords(getGridCoords(c.array() - half_len));
+    std::array<std::size_t, 3> max_coords(getGridCoords(c.array() + half_len));
 
-    tmp_pnt[0] = center[0] + half_len;
-    tmp_pnt[1] = center[1] + half_len;
-    tmp_pnt[2] = center[2] + half_len;
-    std::array<std::size_t, 3> max_coords(getGridCoords(tmp_pnt));
+    std::vector<std::vector<POINT*> const*> pnts;
+    pnts.reserve(
+        (Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(max_coords.data()) -
+         Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(min_coords.data()))
+            .prod());
 
-    std::size_t coords[3], steps0_x_steps1(_n_steps[0] * _n_steps[1]);
-    for (coords[0] = min_coords[0]; coords[0] < max_coords[0] + 1; coords[0]++)
+    std::size_t const steps0_x_steps1(_n_steps[0] * _n_steps[1]);
+    for (std::size_t c0 = min_coords[0]; c0 < max_coords[0] + 1; c0++)
     {
-        for (coords[1] = min_coords[1]; coords[1] < max_coords[1] + 1;
-             coords[1]++)
+        for (std::size_t c1 = min_coords[1]; c1 < max_coords[1] + 1; c1++)
         {
-            const std::size_t coords0_p_coords1_x_steps0(
-                coords[0] + coords[1] * _n_steps[0]);
-            for (coords[2] = min_coords[2]; coords[2] < max_coords[2] + 1;
-                 coords[2]++)
+            const std::size_t coords0_p_coords1_x_steps0(c0 + c1 * _n_steps[0]);
+            for (std::size_t c2 = min_coords[2]; c2 < max_coords[2] + 1; c2++)
             {
                 pnts.push_back(
                     &(_grid_cell_nodes_map[coords0_p_coords1_x_steps0 +
-                                           coords[2] * steps0_x_steps1]));
+                                           c2 * steps0_x_steps1]));
             }
         }
     }
@@ -284,31 +280,34 @@ Grid<POINT>::getPntVecsOfGridCellsIntersectingCube(P const& center,
 }
 
 template <typename POINT>
-void Grid<POINT>::getPntVecsOfGridCellsIntersectingCuboid(
-    Eigen::Vector3d const& min_pnt,
-    Eigen::Vector3d const& max_pnt,
-    std::vector<std::vector<POINT*> const*>& pnts) const
+std::vector<std::vector<POINT*> const*>
+Grid<POINT>::getPntVecsOfGridCellsIntersectingCuboid(
+    Eigen::Vector3d const& min_pnt, Eigen::Vector3d const& max_pnt) const
 {
     std::array<std::size_t, 3> min_coords(getGridCoords(min_pnt));
     std::array<std::size_t, 3> max_coords(getGridCoords(max_pnt));
 
-    std::size_t coords[3], steps0_x_steps1(_n_steps[0] * _n_steps[1]);
-    for (coords[0] = min_coords[0]; coords[0] < max_coords[0] + 1; coords[0]++)
+    std::vector<std::vector<POINT*> const*> pnts;
+    pnts.reserve(
+        (Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(max_coords.data()) -
+         Eigen::Map<Eigen::Matrix<std::size_t, 3, 1>>(min_coords.data()))
+            .prod());
+
+    std::size_t const steps0_x_steps1(_n_steps[0] * _n_steps[1]);
+    for (std::size_t c0 = min_coords[0]; c0 < max_coords[0] + 1; c0++)
     {
-        for (coords[1] = min_coords[1]; coords[1] < max_coords[1] + 1;
-             coords[1]++)
+        for (std::size_t c1 = min_coords[1]; c1 < max_coords[1] + 1; c1++)
         {
-            const std::size_t coords0_p_coords1_x_steps0(
-                coords[0] + coords[1] * _n_steps[0]);
-            for (coords[2] = min_coords[2]; coords[2] < max_coords[2] + 1;
-                 coords[2]++)
+            const std::size_t coords0_p_coords1_x_steps0(c0 + c1 * _n_steps[0]);
+            for (std::size_t c2 = min_coords[2]; c2 < max_coords[2] + 1; c2++)
             {
                 pnts.push_back(
                     &(_grid_cell_nodes_map[coords0_p_coords1_x_steps0 +
-                                           coords[2] * steps0_x_steps1]));
+                                           c2 * steps0_x_steps1]));
             }
         }
     }
+    return pnts;
 }
 
 #ifndef NDEBUG
@@ -421,27 +420,21 @@ std::array<std::size_t, 3> Grid<POINT>::getGridCoords(T const& pnt) const
 {
     auto const& min_point{getMinPoint()};
     auto const& max_point{getMinPoint()};
-    std::array<std::size_t, 3> coords{};
+    std::array<std::size_t, 3> coords{0, 0, 0};
     for (std::size_t k(0); k < 3; k++)
     {
         if (pnt[k] < min_point[k])
         {
-            coords[k] = 0;
+            continue;
         }
-        else
+        if (pnt[k] >= max_point[k])
         {
-            if (pnt[k] >= max_point[k])
-            {
-                coords[k] = _n_steps[k] - 1;
-            }
-            else
-            {
-                coords[k] = static_cast<std::size_t>(
-                    std::floor((pnt[k] - min_point[k])) /
-                    std::nextafter(_step_sizes[k],
-                                   std::numeric_limits<double>::max()));
-            }
+            coords[k] = _n_steps[k] - 1;
+            continue;
         }
+        coords[k] = static_cast<std::size_t>(
+            std::floor((pnt[k] - min_point[k])) /
+            std::nextafter(_step_sizes[k], std::numeric_limits<double>::max()));
     }
     return coords;
 }
diff --git a/GeoLib/PointVec.cpp b/GeoLib/PointVec.cpp
index 15b428ba699cb6c824b0599ac51f938465d74efc..2580e03cb860a90ac07f65c24ad1aa4f57773291 100644
--- a/GeoLib/PointVec.cpp
+++ b/GeoLib/PointVec.cpp
@@ -1,7 +1,5 @@
 /**
  * \file
- * \author Thomas Fischeror
- * \date   2010-06-11
  * \brief  Implementation of the PointVec class.
  *
  * \copyright
diff --git a/GeoLib/PointVec.h b/GeoLib/PointVec.h
index ebea781a1b6a271b50714c9f0674eff60417674f..e42c03d7623dde949f19de095d722e4517bd0f21 100644
--- a/GeoLib/PointVec.h
+++ b/GeoLib/PointVec.h
@@ -22,7 +22,6 @@
 #include "AABB.h"
 #include "OctTree.h"
 #include "Point.h"
-#include "Station.h"
 #include "TemplateVec.h"
 
 namespace GeoLib
diff --git a/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp b/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
index 5c31760f63b915e15b2350ad00bd782aea36c921..3504a8efec15fbd15c1cd8f71e3987c2e713515f 100644
--- a/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
+++ b/MeshLib/MeshEditing/Mesh2MeshPropertyInterpolation.cpp
@@ -116,9 +116,9 @@ void Mesh2MeshPropertyInterpolation::interpolatePropertiesForMesh(
             dest_element.getNodes() + dest_element.getNumberOfBaseNodes());
 
         // request "interesting" nodes from grid
-        std::vector<std::vector<MeshLib::Node*> const*> nodes;
-        src_grid.getPntVecsOfGridCellsIntersectingCuboid(
-            elem_aabb.getMinPoint(), elem_aabb.getMaxPoint(), nodes);
+        std::vector<std::vector<MeshLib::Node*> const*> const nodes =
+            src_grid.getPntVecsOfGridCellsIntersectingCuboid(
+                elem_aabb.getMinPoint(), elem_aabb.getMaxPoint());
 
         std::size_t cnt(0);
         double average_value(0.0);
diff --git a/Tests/GeoLib/TestPointToStationConversion.cpp b/Tests/GeoLib/TestPointToStationConversion.cpp
index 4be5acb588e9f9e98c667970280ff27089b46596..c2263e1f00837079320b9038cac2c12853521a8f 100644
--- a/Tests/GeoLib/TestPointToStationConversion.cpp
+++ b/Tests/GeoLib/TestPointToStationConversion.cpp
@@ -1,6 +1,6 @@
 /**
  * \file
- * \brief Tests for GeoLib::GEOOnjects::geoPointsToStations()
+ * \brief Tests for geoPointsToStations()
  *
  * \copyright
  * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
@@ -33,7 +33,7 @@ TEST(GeoLib, PointToStationConversion)
     // converting only unused points
     std::string stn_orphaned_pnts("Orphaned Points");
     int const res_orphaned_pnts =
-        geo_obj.geoPointsToStations(geo_name, stn_orphaned_pnts, true);
+        geoPointsToStations(geo_obj, geo_name, stn_orphaned_pnts, true);
     EXPECT_EQ(0, res_orphaned_pnts);
     auto const* const stns = geo_obj.getStationVec(stn_orphaned_pnts);
     assert(stns != nullptr);
@@ -49,7 +49,7 @@ TEST(GeoLib, PointToStationConversion)
     // converting all points
     std::string stn_all_pnts("All Points");
     int const res_all_pnts =
-        geo_obj.geoPointsToStations(geo_name, stn_all_pnts, false);
+        geoPointsToStations(geo_obj, geo_name, stn_all_pnts, false);
     EXPECT_EQ(0, res_all_pnts);
     EXPECT_EQ(exp_all_pnts, geo_obj.getStationVec(stn_all_pnts)->size());
 }