diff --git a/BaseLib/Subdivision.cpp b/BaseLib/Subdivision.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a2141d04d7728022f8e9ec9eacd5f343959a9e38
--- /dev/null
+++ b/BaseLib/Subdivision.cpp
@@ -0,0 +1,70 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "Subdivision.h"
+
+#include <BaseLib/Error.h>
+
+namespace BaseLib
+{
+GradualSubdivision::GradualSubdivision(const double L,
+                                       const double dL0,
+                                       const double max_dL,
+                                       const double multiplier)
+    : _length(L), _dL0(dL0), _max_dL(max_dL), _multiplier(multiplier)
+{
+    // Check if accumulated subdivisions can ever sum up to length.
+    // Cf. geometric series formula.
+    if (multiplier < 1.0 && dL0 / (1.0 - multiplier) < L) {
+        OGS_FATAL(
+            "Using dL0=%g and multiplier=%g the generated subdivisions can not "
+            "sum up to a total length of %g.",
+            dL0,
+            multiplier,
+            L);
+    }
+}
+
+GradualSubdivisionFixedNum::GradualSubdivisionFixedNum(
+    const double L, const std::size_t num_subdivisions, const double multiplier)
+    : _length{L}, _num_subdivisions{num_subdivisions}, _multiplier{multiplier}
+{
+}
+
+std::vector<double> GradualSubdivisionFixedNum::operator()() const
+{
+    std::vector<double> subdivisions;
+    subdivisions.reserve(_num_subdivisions + 1);
+    subdivisions.push_back(0.0);
+    auto const q = _multiplier;
+
+    if (q == 1.0) {
+        double const dx = _length / _num_subdivisions;
+
+        for (std::size_t i = 1; i < _num_subdivisions; ++i) {
+            subdivisions.push_back(dx * i);
+        }
+    } else {
+        // compute initial subdivision size
+        auto const a =
+            _length * (q - 1.0) / (std::pow(q, _num_subdivisions) - 1.0);
+
+        double qi = q;  // q^i
+        for (std::size_t i = 1; i < _num_subdivisions; ++i) {
+            subdivisions.push_back(a * (qi - 1.0) / (q - 1.0));
+            qi *= q;
+        }
+    }
+
+    subdivisions.push_back(_length);
+
+    return subdivisions;
+}
+
+}  // namespace BaseLib
diff --git a/BaseLib/Subdivision.h b/BaseLib/Subdivision.h
index 625c1efe76800ea712466366a799cd71dfa46db1..f6c87e869849466482acf4d33c3f1cea957e3ae5 100644
--- a/BaseLib/Subdivision.h
+++ b/BaseLib/Subdivision.h
@@ -10,6 +10,7 @@
 #pragma once
 
 #include <cmath>
+#include <vector>
 
 namespace BaseLib
 {
@@ -73,8 +74,7 @@ public:
             const double L,
             const double dL0,
             const double max_dL,
-            const double multiplier)
-    : _length(L), _dL0(dL0), _max_dL(max_dL), _multiplier(multiplier) {}
+            const double multiplier);
 
     /// Returns a vector of subdivided points
     std::vector<double> operator()() const override
@@ -106,4 +106,31 @@ private:
     const double _multiplier;
 };
 
+/**
+ * Gradual subdivision operator with a constant _multiplier.
+ *
+ * In this class the number of subdivisions is known a priori.
+ */
+class GradualSubdivisionFixedNum : public ISubdivision
+{
+public:
+    /**
+     * Constructor
+     * @param L                 total length to be subdivided
+     * @param num_subdivisions  number of subdivisions to generate
+     * @param multiplier        multiplier to cell length
+     */
+    GradualSubdivisionFixedNum(const double L,
+                               const std::size_t num_subdivisions,
+                               const double multiplier);
+
+    /// Returns a vector of subdivided points
+    std::vector<double> operator()() const override;
+
+private:
+    const double _length;
+    const std::size_t _num_subdivisions;
+    const double _multiplier;
+};
+
 } // BaseLib