diff --git a/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp b/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp
index 04ef2e9a5f8276fa7c4bd79c0e9a4b528fc0dd05..06f65f8cc31942acc0af052053126a9bc07cfa9e 100644
--- a/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp
+++ b/Applications/Utils/ModelPreparation/PartitionMesh/PartitionMesh.cpp
@@ -133,7 +133,7 @@ int main(int argc, char* argv[])
                 std::ofstream os(name);
                 if (!os)
                     OGS_FATAL("Couldn't open file '%s' for writing.",
-                              name.c_str())
+                              name.c_str());
                 for (std::size_t n(0); n < number; ++n)
                     os << "0\n";
             };
diff --git a/Applications/Utils/SimpleMeshCreation/generateStructuredMesh.cpp b/Applications/Utils/SimpleMeshCreation/generateStructuredMesh.cpp
index cdc7ce9fbd99be84da828bf881a844d93771a8e3..7d26ec3c325644bb65ac284c17bb79e76dfd9197 100644
--- a/Applications/Utils/SimpleMeshCreation/generateStructuredMesh.cpp
+++ b/Applications/Utils/SimpleMeshCreation/generateStructuredMesh.cpp
@@ -16,6 +16,7 @@
 #include "Applications/ApplicationsLib/LogogSetup.h"
 
 #include "BaseLib/BuildInfo.h"
+#include "BaseLib/Error.h"
 #include "BaseLib/Subdivision.h"
 #include "BaseLib/TCLAPCustomOutput.h"
 
@@ -179,17 +180,27 @@ int main (int argc, char* argv[])
     vec_div.reserve(dim);
     for (unsigned i=0; i<dim; i++)
     {
-        if (vec_ndivArg[i]->isSet())
-        {
+        if (vec_multiArg[i]->isSet()) {
+            if (vec_ndivArg[i]->isSet()) {
+                // number of partitions in direction is specified
+                if (vec_d0Arg[i]->isSet()) {
+                    OGS_FATAL(
+                        "Specifying all of --m?, --d?0 and --n? for coordinate "
+                        "\"?\" is not supported.");
+                }
+                vec_div.emplace_back(new BaseLib::GradualSubdivisionFixedNum(
+                    length[i], vec_ndivArg[i]->getValue(),
+                    vec_multiArg[i]->getValue()));
+
+            } else {
+                vec_div.emplace_back(new BaseLib::GradualSubdivision(
+                    length[i], vec_d0Arg[i]->getValue(),
+                    vec_dMaxArg[i]->getValue(), vec_multiArg[i]->getValue()));
+            }
+        } else {
             vec_div.emplace_back(
                 new BaseLib::UniformSubdivision(length[i], n_subdivision[i]));
         }
-        else
-        {
-            vec_div.emplace_back(new BaseLib::GradualSubdivision(
-                length[i], vec_d0Arg[i]->getValue(), vec_dMaxArg[i]->getValue(),
-                vec_multiArg[i]->getValue()));
-        }
     }
 
     // generate a mesh
diff --git a/BaseLib/Error.h b/BaseLib/Error.h
index 517c051bca2218826b779669fe4f7c7fc67403af..fd5503c9d4240a833f61bf85163e7dd4873d3284 100644
--- a/BaseLib/Error.h
+++ b/BaseLib/Error.h
@@ -13,25 +13,51 @@
 
 #include <cstdlib>
 #include <logog/include/logog.hpp>
-#include "StringTools.h"
 
-#define OGS_FATAL(fmt, ...)\
-    {\
-    ERR("%s", BaseLib::format(fmt, ##__VA_ARGS__).data());\
-    std::abort();\
-    }
+namespace BaseLib
+{
+namespace detail
+{
+template <typename Msg>
+[[noreturn]] bool error_impl(Msg&& msg)
+{
+    ERR("%s", msg.data());
+    std::abort();
+}
+
+}  // namespace detail
+
+}  // namespace BaseLib
 
-#else // OGS_FATAL_ABORT
+#else  // OGS_FATAL_ABORT
 
 #include <stdexcept>
+
+namespace BaseLib
+{
+namespace detail
+{
+template <typename Msg>
+[[noreturn]] bool error_impl(Msg&& msg)
+{
+    throw std::runtime_error(std::forward<Msg>(msg));
+}
+
+}  // namespace detail
+
+}  // namespace BaseLib
+
+#endif  // OGS_FATAL_ABORT
+
 #include "FileTools.h"
 #include "StringTools.h"
 
 #define OGS_STR(x) #x
 #define OGS_STRINGIFY(x) OGS_STR(x)
-#define OGS_LOCATION " at " + BaseLib::extractBaseName(__FILE__) + ", line " OGS_STRINGIFY(__LINE__)
-
-#define OGS_FATAL(fmt, ...)\
-    throw std::runtime_error(BaseLib::format(fmt, ##__VA_ARGS__) + OGS_LOCATION);
+#define OGS_LOCATION                              \
+    " at " + BaseLib::extractBaseName(__FILE__) + \
+        ", line " OGS_STRINGIFY(__LINE__)
 
-#endif // OGS_FATAL_ABORT
+#define OGS_FATAL(fmt, ...)                                           \
+    BaseLib::detail::error_impl(BaseLib::format(fmt, ##__VA_ARGS__) + \
+                                OGS_LOCATION)
diff --git a/BaseLib/Subdivision.cpp b/BaseLib/Subdivision.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b9ba827c5bc01619bf6da6eff0be2f6149588d70
--- /dev/null
+++ b/BaseLib/Subdivision.cpp
@@ -0,0 +1,96 @@
+/**
+ * \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 <algorithm>
+#include <cmath>
+
+#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);
+    }
+}
+
+std::vector<double> GradualSubdivision::operator()() const
+{
+    std::vector<double> vec_x;
+
+    double x = 0;
+    unsigned i = 0;
+    do {
+        vec_x.push_back(x);
+        x += std::min(_max_dL,
+                      _dL0 * std::pow(_multiplier, static_cast<double>(i)));
+        i++;
+    } while (x < _length);
+
+    if (vec_x.back() < _length) {
+        double last_dx = vec_x[vec_x.size() - 1] - vec_x[vec_x.size() - 2];
+        if (_length - vec_x.back() < last_dx)
+            vec_x[vec_x.size() - 1] = _length;
+        else
+            vec_x.push_back(_length);
+    }
+    return vec_x;
+}
+
+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..442fb514af91b15994e98aa057a3e4a41aaab0c6 100644
--- a/BaseLib/Subdivision.h
+++ b/BaseLib/Subdivision.h
@@ -9,7 +9,7 @@
 
 #pragma once
 
-#include <cmath>
+#include <vector>
 
 namespace BaseLib
 {
@@ -73,31 +73,10 @@ 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
-    {
-        std::vector<double> vec_x;
-
-        double x = 0;
-        unsigned i=0;
-        do {
-            vec_x.push_back(x);
-            x += std::min(_max_dL, _dL0*std::pow(_multiplier, static_cast<double>(i)));
-            i++;
-        } while (x<_length);
-
-        if (vec_x.back() < _length) {
-            double last_dx = vec_x[vec_x.size()-1] - vec_x[vec_x.size()-2];
-            if (_length-vec_x.back()<last_dx)
-                vec_x[vec_x.size()-1] = _length;
-            else
-                vec_x.push_back(_length);
-        }
-        return vec_x;
-    }
+    std::vector<double> operator()() const override;
 
 private:
     const double _length;
@@ -106,4 +85,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
diff --git a/MaterialLib/Adsorption/ReactionInert.h b/MaterialLib/Adsorption/ReactionInert.h
index 3e3a05052c1a970f99b3d03ee18b793342b0d5a2..f6cd26302fb5a4ed9e8455e57c1cfdb8257f8ec6 100644
--- a/MaterialLib/Adsorption/ReactionInert.h
+++ b/MaterialLib/Adsorption/ReactionInert.h
@@ -28,7 +28,7 @@ public:
     double getReactionRate(const double /*p_Ads*/, const double /*T_Ads*/, const double /*M_Ads*/,
                              const double /*loading*/) const override
     {
-        OGS_FATAL("Method getReactionRate() should never be called directly")
+        OGS_FATAL("Method getReactionRate() should never be called directly");
     }
 };
 
diff --git a/MaterialLib/PorousMedium/Permeability/createPermeabilityModel.cpp b/MaterialLib/PorousMedium/Permeability/createPermeabilityModel.cpp
index 8ffe381b6b3aaabfcde23948b14757dd4bf7dd09..3927950dd2bebd1f454aa6c3ce3530c97ad14dab 100644
--- a/MaterialLib/PorousMedium/Permeability/createPermeabilityModel.cpp
+++ b/MaterialLib/PorousMedium/Permeability/createPermeabilityModel.cpp
@@ -46,7 +46,7 @@ Eigen::MatrixXd createPermeabilityModel(BaseLib::ConfigTree const& config)
         default:
         {
             OGS_FATAL(
-                "Number of values for permeability tensor must be 1, 4 or 9.")
+                "Number of values for permeability tensor must be 1, 4 or 9.");
         }
     }