diff --git a/Applications/FileIO/Gmsh/GMSHInterface.cpp b/Applications/FileIO/Gmsh/GMSHInterface.cpp
index 5f02a6f17a2dcbb778a4688040a9a5002e2eeab7..06ec29143a457db188bd4159385c5dfe22278934 100644
--- a/Applications/FileIO/Gmsh/GMSHInterface.cpp
+++ b/Applications/FileIO/Gmsh/GMSHInterface.cpp
@@ -33,12 +33,13 @@ namespace FileIO
 {
 namespace GMSH
 {
-GMSHInterface::GMSHInterface(GeoLib::GEOObjects& geo_objs,
-                             bool /*include_stations_as_constraints*/,
-                             GMSH::MeshDensityAlgorithm mesh_density_algorithm,
-                             double param1, double param2, std::size_t param3,
-                             std::vector<std::string>& selected_geometries,
-                             bool rotate, bool keep_preprocessed_geometry)
+GMSHInterface::GMSHInterface(
+    GeoLib::GEOObjects& geo_objs, bool /*include_stations_as_constraints*/,
+    GMSH::MeshDensityAlgorithm const mesh_density_algorithm,
+    double const pnt_density, double const station_density,
+    std::size_t const max_pnts_per_leaf,
+    std::vector<std::string> const& selected_geometries, bool const rotate,
+    bool const keep_preprocessed_geometry)
     : _n_lines(0),
       _n_plane_sfc(0),
       _geo_objs(geo_objs),
@@ -48,10 +49,13 @@ GMSHInterface::GMSHInterface(GeoLib::GEOObjects& geo_objs,
 {
     switch (mesh_density_algorithm) {
     case GMSH::MeshDensityAlgorithm::FixedMeshDensity:
-        _mesh_density_strategy = new GMSH::GMSHFixedMeshDensity(param1);
+        _mesh_density_strategy =
+            std::make_unique<GMSH::GMSHFixedMeshDensity>(pnt_density);
         break;
     case GMSH::MeshDensityAlgorithm::AdaptiveMeshDensity:
-        _mesh_density_strategy = new GMSH::GMSHAdaptiveMeshDensity(param1, param2, param3);
+        _mesh_density_strategy =
+            std::make_unique<GMSH::GMSHAdaptiveMeshDensity>(
+                pnt_density, station_density, max_pnts_per_leaf);
         break;
     }
 }
@@ -60,7 +64,6 @@ GMSHInterface::~GMSHInterface()
 {
     for (auto * gmsh_pnt : _gmsh_pnts)
         delete gmsh_pnt;
-    delete _mesh_density_strategy;
     for (auto * polygon_tree : _polygon_tree_list)
         delete polygon_tree;
 }
@@ -116,7 +119,7 @@ int GMSHInterface::writeGMSHInputFile(std::ostream& out)
 
     std::vector<GeoLib::Polyline*> const* merged_plys(
         _geo_objs.getPolylineVec(_gmsh_geo_name));
-    DBUG("GMSHInterface::writeGMSHInputFile(): \t ok.");
+    DBUG("GMSHInterface::writeGMSHInputFile(): Obtained data.");
 
     if (!merged_plys) {
         ERR("GMSHInterface::writeGMSHInputFile(): Did not find any polylines.");
@@ -136,7 +139,7 @@ int GMSHInterface::writeGMSHInputFile(std::ostream& out)
         }
         _polygon_tree_list.push_back(new GMSH::GMSHPolygonTree(
             new GeoLib::PolygonWithSegmentMarker(*polyline), nullptr, _geo_objs,
-            _gmsh_geo_name, _mesh_density_strategy));
+            _gmsh_geo_name, *_mesh_density_strategy));
     }
     DBUG(
         "GMSHInterface::writeGMSHInputFile(): Computed topological hierarchy - "
diff --git a/Applications/FileIO/Gmsh/GMSHInterface.h b/Applications/FileIO/Gmsh/GMSHInterface.h
index 5eecacd72b2ed8f0048591838c566129978b8a01..ae8d68f3265fcd425cffc92640463b4e8d32f02b 100644
--- a/Applications/FileIO/Gmsh/GMSHInterface.h
+++ b/Applications/FileIO/Gmsh/GMSHInterface.h
@@ -52,9 +52,9 @@ public:
      * as constraints
      * @param mesh_density_algorithm one of the mesh density algorithms (\@see
      * enum MeshDensityAlgorithm)
-     * @param param1 parameter that can be used for the mesh density algorithm
-     * @param param2 parameter that can be used for the mesh density algorithm
-     * @param param3 parameter that can be used for the mesh density algorithm
+     * @param pnt_density parameter of the mesh density algorithm
+     * @param station_density parameter of the mesh density algorithm
+     * @param max_pnts_per_leaf parameter of the mesh density algorithm
      * @param selected_geometries vector of names of geometries, that should be
      * employed for mesh generation.
      * @param rotate if the value of the parameter is true then the input points
@@ -66,9 +66,10 @@ public:
     GMSHInterface(GeoLib::GEOObjects& geo_objs,
                   bool include_stations_as_constraints,
                   GMSH::MeshDensityAlgorithm mesh_density_algorithm,
-                  double param1, double param2, std::size_t param3,
-                  std::vector<std::string>& selected_geometries, bool rotate,
-                  bool keep_preprocessed_geometry);
+                  double pnt_density, double station_density,
+                  std::size_t max_pnts_per_leaf,
+                  std::vector<std::string> const& selected_geometries,
+                  bool rotate, bool keep_preprocessed_geometry);
 
     GMSHInterface(GMSHInterface const&) = delete;
     GMSHInterface(GMSHInterface &&) = delete;
@@ -97,20 +98,20 @@ private:
     std::size_t _n_plane_sfc;
 
     GeoLib::GEOObjects & _geo_objs;
-    std::vector<std::string>& _selected_geometries;
+    std::vector<std::string> const& _selected_geometries;
     std::string _gmsh_geo_name;
     std::list<GMSH::GMSHPolygonTree*> _polygon_tree_list;
 
     std::vector<GMSH::GMSHPoint*> _gmsh_pnts;
 
-    GMSH::GMSHMeshDensityStrategy *_mesh_density_strategy;
+    std::unique_ptr<GMSH::GMSHMeshDensityStrategy> _mesh_density_strategy;
     /// Holds the inverse rotation matrix. The matrix is used in writePoints() to
     /// revert the rotation done in writeGMSHInputFile().
     MathLib::DenseMatrix<double> _inverse_rot_mat =
         MathLib::DenseMatrix<double>(3, 3, 0);
     /// Signals if the input points should be rotated or projected to the
     /// \f$x\f$-\f$y\f$-plane
-    bool _rotate = false;
+    bool const _rotate = false;
     bool _keep_preprocessed_geometry = true;
 };
 } // end namespace GMSH
diff --git a/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp b/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
index 7580985583475e1bf931f843fb44a535c938fdb0..e61cdf3ccf1ac8dc95eb0f2b54c4b4493d1d572a 100644
--- a/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
+++ b/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
@@ -24,13 +24,15 @@ namespace FileIO
 {
 namespace GMSH
 {
-
 GMSHPolygonTree::GMSHPolygonTree(GeoLib::PolygonWithSegmentMarker* polygon,
-                GMSHPolygonTree* parent,
-                GeoLib::GEOObjects &geo_objs, std::string const& geo_name,
-                GMSHMeshDensityStrategy * mesh_density_strategy) :
-    GeoLib::SimplePolygonTree(polygon, parent), _geo_objs(geo_objs), _geo_name(geo_name),
-    _mesh_density_strategy(mesh_density_strategy)
+                                 GMSHPolygonTree* parent,
+                                 GeoLib::GEOObjects& geo_objs,
+                                 std::string const& geo_name,
+                                 GMSHMeshDensityStrategy& mesh_density_strategy)
+    : GeoLib::SimplePolygonTree(polygon, parent),
+      _geo_objs(geo_objs),
+      _geo_name(geo_name),
+      _mesh_density_strategy(mesh_density_strategy)
 {}
 
 GMSHPolygonTree::~GMSHPolygonTree()
@@ -222,7 +224,9 @@ void GMSHPolygonTree::checkIntersectionsSegmentExistingPolylines(
 
 void GMSHPolygonTree::initMeshDensityStrategy()
 {
-    if (dynamic_cast<GMSHAdaptiveMeshDensity*> (_mesh_density_strategy)) {
+    if (auto* adaptive_mesh_density =
+            dynamic_cast<GMSHAdaptiveMeshDensity*>(&_mesh_density_strategy))
+    {
         // collect points
         std::vector<GeoLib::Point const*> pnts;
         const std::size_t n_pnts_polygon (_node_polygon->getNumberOfPoints());
@@ -240,12 +244,12 @@ void GMSHPolygonTree::initMeshDensityStrategy()
         }
 
         // give collected points to the mesh density strategy
-        _mesh_density_strategy->initialize(pnts);
+        adaptive_mesh_density->initialize(pnts);
         // insert constraints
-        dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)->addPoints(_stations);
+        adaptive_mesh_density->addPoints(_stations);
         std::vector<GeoLib::Point const*> stations;
         getStationsInsideSubPolygons(stations);
-        dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)->addPoints(stations);
+        adaptive_mesh_density->addPoints(stations);
     }
 }
 
@@ -259,7 +263,7 @@ void GMSHPolygonTree::createGMSHPoints(std::vector<GMSHPoint*> & gmsh_pnts) cons
         if (gmsh_pnts[id] != nullptr)
             continue;
         gmsh_pnts[id] = new GMSHPoint(
-            *pnt, id, _mesh_density_strategy->getMeshDensityAtPoint(pnt));
+            *pnt, id, _mesh_density_strategy.getMeshDensityAtPoint(pnt));
     }
 
     const std::size_t n_plys(_plys.size());
@@ -274,7 +278,7 @@ void GMSHPolygonTree::createGMSHPoints(std::vector<GMSHPoint*> & gmsh_pnts) cons
                 GeoLib::Point const*const pnt(_plys[k]->getPoint(j));
                 gmsh_pnts[id] = new GMSHPoint(
                     *pnt, id,
-                    _mesh_density_strategy->getMeshDensityAtPoint(pnt));
+                    _mesh_density_strategy.getMeshDensityAtPoint(pnt));
             }
         }
     }
@@ -359,7 +363,7 @@ void GMSHPolygonTree::writeStations(std::size_t & pnt_id_offset, std::size_t sfc
     for (auto const* station : _stations) {
         out << "Point(" << pnt_id_offset << ") = {" << (*station)[0] << ", "
             << (*station)[1] << ", 0.0, "
-            << _mesh_density_strategy->getMeshDensityAtStation(station)
+            << _mesh_density_strategy.getMeshDensityAtStation(station)
             << "}; // Station "
             << static_cast<GeoLib::Station const*>(station)->getName() << " \n";
         out << "Point { " << pnt_id_offset << " } In Surface { " << sfc_number << " };\n";
@@ -369,14 +373,18 @@ void GMSHPolygonTree::writeStations(std::size_t & pnt_id_offset, std::size_t sfc
 
 void GMSHPolygonTree::writeAdditionalPointData(std::size_t & pnt_id_offset, std::size_t sfc_number, std::ostream& out) const
 {
-    if (dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)) {
+    if (auto* adaptive_mesh_density =
+            dynamic_cast<GMSHAdaptiveMeshDensity*>(&_mesh_density_strategy))
+    {
         std::vector<GeoLib::Point*> steiner_pnts;
-        dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)->getSteinerPoints(steiner_pnts, 0);
+        adaptive_mesh_density->getSteinerPoints(steiner_pnts, 0);
         const std::size_t n(steiner_pnts.size());
         for (std::size_t k(0); k<n; k++) {
             if (_node_polygon->isPntInPolygon(*(steiner_pnts[k]))) {
                 out << "Point(" << pnt_id_offset + k << ") = {" << (*(steiner_pnts[k]))[0] << "," << (*(steiner_pnts[k]))[1] << ", 0.0, ";
-                out << _mesh_density_strategy->getMeshDensityAtPoint(steiner_pnts[k]) << "};\n";
+                out << _mesh_density_strategy.getMeshDensityAtPoint(
+                           steiner_pnts[k])
+                    << "};\n";
                 out << "Point { " << pnt_id_offset + k << " } In Surface { " << sfc_number << " };\n";
             }
             delete steiner_pnts[k];
@@ -385,10 +393,12 @@ void GMSHPolygonTree::writeAdditionalPointData(std::size_t & pnt_id_offset, std:
     }
 
 #ifndef NDEBUG
-    if (dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)) {
+    if (auto* adaptive_mesh_density =
+            dynamic_cast<GMSHAdaptiveMeshDensity*>(&_mesh_density_strategy))
+    {
         auto pnts = std::make_unique<std::vector<GeoLib::Point*>>();
         auto plys = std::make_unique<std::vector<GeoLib::Polyline*>>();
-        dynamic_cast<GMSHAdaptiveMeshDensity*>(_mesh_density_strategy)->getQuadTreeGeometry(*pnts, *plys);
+        adaptive_mesh_density->getQuadTreeGeometry(*pnts, *plys);
         std::string quad_tree_geo("QuadTree");
         _geo_objs.addPointVec(std::move(pnts), quad_tree_geo);
         std::vector<std::size_t> const& id_map ((_geo_objs.getPointVecObj(quad_tree_geo))->getIDMap());
diff --git a/Applications/FileIO/Gmsh/GMSHPolygonTree.h b/Applications/FileIO/Gmsh/GMSHPolygonTree.h
index dac4afa61abc6c44be4aec740f1091d3108f3724..257900c3cf3630e13db121b361a696be3232fe37 100644
--- a/Applications/FileIO/Gmsh/GMSHPolygonTree.h
+++ b/Applications/FileIO/Gmsh/GMSHPolygonTree.h
@@ -34,9 +34,10 @@ namespace GMSH
 
 class GMSHPolygonTree: public GeoLib::SimplePolygonTree {
 public:
-    GMSHPolygonTree(GeoLib::PolygonWithSegmentMarker* polygon, GMSHPolygonTree * parent,
-                    GeoLib::GEOObjects &geo_objs, std::string const& geo_name,
-                    GMSHMeshDensityStrategy * mesh_density_strategy);
+    GMSHPolygonTree(GeoLib::PolygonWithSegmentMarker* polygon,
+                    GMSHPolygonTree* parent, GeoLib::GEOObjects& geo_objs,
+                    std::string const& geo_name,
+                    GMSHMeshDensityStrategy& mesh_density_strategy);
     ~GMSHPolygonTree() override;
 
     /** Mark the segments shared by several polygons. */
@@ -96,7 +97,7 @@ private:
     std::vector<GeoLib::PolylineWithSegmentMarker*> _plys;
     std::vector<GMSHLine*> _gmsh_lines_for_constraints;
 
-    GMSHMeshDensityStrategy * _mesh_density_strategy;
+    GMSHMeshDensityStrategy& _mesh_density_strategy;
 };
 
 }  // end namespace GMSH
diff --git a/GeoLib/PointVec.h b/GeoLib/PointVec.h
index 1aad6c17f5d64918538cb8dd264c95c965cc9173..ad936eef806a76ba34296febe90988d4d78c4bf8 100644
--- a/GeoLib/PointVec.h
+++ b/GeoLib/PointVec.h
@@ -127,6 +127,17 @@ private:
     // this way the compiler does not create a (possible unwanted) assignment operator
     PointVec& operator= (const PointVec& rhs);
 
+	/**
+	 * Inserts the instance of the Point into internal data structures
+	 * (@see TemplateVec::_data_vec, _pnt_id_map) if and only if there
+	 * does not exist a point with the same coordinates (upto
+	 * std::numeric_limits<double>::epsilon()). In case there exists
+	 * already a point with the same coordinates the given pnt-object
+	 * will be deleted!
+	 * @param pnt  Pointer to GeooLib::Point instance
+	 * @return either the new id or the id of the existing point with
+	 * the same coordinates
+	 */
     std::size_t uniqueInsert (Point* pnt);
 
     /** the type of the point (\sa enum PointType) */