diff --git a/NumLib/DOF/ComputeSparsityPattern.cpp b/NumLib/DOF/ComputeSparsityPattern.cpp
index 4bcbcb62562fe596652cf90da833a8fb81529477..166ecc0efe8654b430e0df43e723ec777ea0819d 100644
--- a/NumLib/DOF/ComputeSparsityPattern.cpp
+++ b/NumLib/DOF/ComputeSparsityPattern.cpp
@@ -19,8 +19,7 @@
 #include "MeshLib/NodePartitionedMesh.h"
 
 GlobalSparsityPattern computeSparsityPatternPETSc(
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    MeshLib::Mesh const& mesh)
+    NumLib::LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
 {
     assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(&mesh));
     auto const& npmesh =
diff --git a/NumLib/DOF/DOFTableUtil.cpp b/NumLib/DOF/DOFTableUtil.cpp
index 1364d5b9b116b21821f7593776c19cb4ae76df30..7879eb2c5c5f5b43e478d44e2acb40baeb87d85e 100644
--- a/NumLib/DOF/DOFTableUtil.cpp
+++ b/NumLib/DOF/DOFTableUtil.cpp
@@ -9,6 +9,7 @@
  */
 
 #include "DOFTableUtil.h"
+
 #include <cassert>
 
 namespace NumLib
@@ -97,7 +98,7 @@ double getNonGhostNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
                               node_id};
 
     auto const index = dof_table.getGlobalIndex(l, global_component_id);
-    assert (index != NumLib::MeshComponentMap::nop);
+    assert(index != NumLib::MeshComponentMap::nop);
 
     if (index < 0)
     {  // ghost node value
@@ -116,7 +117,7 @@ double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
                               node_id};
 
     auto const index = dof_table.getGlobalIndex(l, global_component_id);
-    assert (index != NumLib::MeshComponentMap::nop);
+    assert(index != NumLib::MeshComponentMap::nop);
 
     return x.get(index);
 }
diff --git a/NumLib/DOF/GlobalMatrixProviders.cpp b/NumLib/DOF/GlobalMatrixProviders.cpp
index 10c8628513530bda39ce8e3ddad39fc3d181d5bc..690ddf851eb03183c3fc0fa9a367adcc92f76ba8 100644
--- a/NumLib/DOF/GlobalMatrixProviders.cpp
+++ b/NumLib/DOF/GlobalMatrixProviders.cpp
@@ -14,26 +14,25 @@
 
 #include "SimpleMatrixVectorProvider.h"
 
-
 // Initializes the static members of the structs in the header file
 // associated with this file.
-#define INITIALIZE_GLOBAL_MATRIX_VECTOR_PROVIDER(VARNAME) \
+#define INITIALIZE_GLOBAL_MATRIX_VECTOR_PROVIDER(VARNAME)               \
     static std::unique_ptr<NumLib::SimpleMatrixVectorProvider> VARNAME{ \
-        new NumLib::SimpleMatrixVectorProvider}; \
-    \
-    namespace NumLib { \
-    VectorProvider& GlobalVectorProvider::provider = *(VARNAME); \
-    \
-    MatrixProvider& GlobalMatrixProvider::provider = *(VARNAME); \
+        new NumLib::SimpleMatrixVectorProvider};                        \
+                                                                        \
+    namespace NumLib                                                    \
+    {                                                                   \
+    VectorProvider& GlobalVectorProvider::provider = *(VARNAME);        \
+                                                                        \
+    MatrixProvider& GlobalMatrixProvider::provider = *(VARNAME);        \
     }
 
 INITIALIZE_GLOBAL_MATRIX_VECTOR_PROVIDER(globalSetupGlobalMatrixVectorProvider)
 
-
 namespace NumLib
 {
 void cleanupGlobalMatrixProviders()
 {
     globalSetupGlobalMatrixVectorProvider.reset();
 }
-}
+}  // namespace NumLib
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.cpp b/NumLib/DOF/LocalToGlobalIndexMap.cpp
index 1ef6b39636184fc61603b78e0b169b886116e129..cd1c61e4d77036e137896affeb936f33c260c4f6 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.cpp
+++ b/NumLib/DOF/LocalToGlobalIndexMap.cpp
@@ -16,16 +16,14 @@
 
 namespace NumLib
 {
-
 namespace
 {
-
 // Make the cumulative sum of an array, which starts with zero
 template <typename T>
 std::vector<T> to_cumulative(std::vector<T> const& vec)
 {
-    std::vector<T> result(vec.size()+1, 0);
-    std::partial_sum(vec.begin(), vec.end(), result.begin()+1);
+    std::vector<T> result(vec.size() + 1, 0);
+    std::partial_sum(vec.begin(), vec.end(), result.begin() + 1);
 
     return result;
 }
@@ -44,7 +42,8 @@ void LocalToGlobalIndexMap::findGlobalIndicesWithElementID(
     std::vector<MeshLib::Node*> const& nodes, std::size_t const mesh_id,
     const int comp_id, const int comp_id_write)
 {
-    std::unordered_set<MeshLib::Node*> const set_nodes(nodes.begin(), nodes.end());
+    std::unordered_set<MeshLib::Node*> const set_nodes(nodes.begin(),
+                                                       nodes.end());
 
     // For each element find the global indices for node/element
     // components.
@@ -54,15 +53,16 @@ void LocalToGlobalIndexMap::findGlobalIndicesWithElementID(
         indices.reserve((*e)->getNumberOfNodes());
 
         for (auto* n = (*e)->getNodes();
-             n < (*e)->getNodes()+(*e)->getNumberOfNodes(); ++n)
+             n < (*e)->getNodes() + (*e)->getNumberOfNodes();
+             ++n)
         {
             // Check if the element's node is in the given list of nodes.
             if (set_nodes.find(*n) == set_nodes.end())
             {
                 continue;
             }
-            MeshLib::Location l(
-                mesh_id, MeshLib::MeshItemType::Node, (*n)->getID());
+            MeshLib::Location l(mesh_id, MeshLib::MeshItemType::Node,
+                                (*n)->getID());
             indices.push_back(_mesh_component_map.getGlobalIndex(l, comp_id));
         }
 
@@ -79,7 +79,8 @@ void LocalToGlobalIndexMap::findGlobalIndices(
 {
     _rows.resize(std::distance(first, last), _mesh_subsets.size());
 
-    std::unordered_set<MeshLib::Node*> const set_nodes(nodes.begin(), nodes.end());
+    std::unordered_set<MeshLib::Node*> const set_nodes(nodes.begin(),
+                                                       nodes.end());
 
     // For each element find the global indices for node/element
     // components.
@@ -90,15 +91,16 @@ void LocalToGlobalIndexMap::findGlobalIndices(
         indices.reserve((*e)->getNumberOfNodes());
 
         for (auto* n = (*e)->getNodes();
-             n < (*e)->getNodes() + (*e)->getNumberOfNodes(); ++n)
+             n < (*e)->getNodes() + (*e)->getNumberOfNodes();
+             ++n)
         {
             // Check if the element's node is in the given list of nodes.
             if (set_nodes.find(*n) == set_nodes.end())
             {
                 continue;
             }
-            MeshLib::Location l(
-                mesh_id, MeshLib::MeshItemType::Node, (*n)->getID());
+            MeshLib::Location l(mesh_id, MeshLib::MeshItemType::Node,
+                                (*n)->getID());
             auto const global_index =
                 _mesh_component_map.getGlobalIndex(l, comp_id);
             if (global_index == std::numeric_limits<GlobalIndexType>::max())
@@ -130,10 +132,12 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
       _variable_component_offsets(to_cumulative(vec_var_n_components))
 {
     // For each element of that MeshSubset save a line of global indices.
-    for (int variable_id = 0; variable_id < static_cast<int>(vec_var_n_components.size());
+    for (int variable_id = 0;
+         variable_id < static_cast<int>(vec_var_n_components.size());
          ++variable_id)
     {
-        for (int component_id = 0; component_id < static_cast<int>(vec_var_n_components[variable_id]);
+        for (int component_id = 0;
+             component_id < static_cast<int>(vec_var_n_components[variable_id]);
              ++component_id)
         {
             auto const global_component_id =
@@ -164,7 +168,7 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
 
     // _rows should be resized based on an element ID
     std::size_t max_elem_id = 0;
-    for (std::vector<MeshLib::Element*>const* elements : vec_var_elements)
+    for (std::vector<MeshLib::Element*> const* elements : vec_var_elements)
     {
         for (auto e : *elements)
         {
@@ -173,11 +177,14 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
     }
     _rows.resize(max_elem_id + 1, _mesh_subsets.size());
 
-    for (int variable_id = 0; variable_id < static_cast<int>(vec_var_n_components.size());
+    for (int variable_id = 0;
+         variable_id < static_cast<int>(vec_var_n_components.size());
          ++variable_id)
     {
-        std::vector<MeshLib::Element*> const& var_elements = *vec_var_elements[variable_id];
-        for (int component_id = 0; component_id < static_cast<int>(vec_var_n_components[variable_id]);
+        std::vector<MeshLib::Element*> const& var_elements =
+            *vec_var_elements[variable_id];
+        for (int component_id = 0;
+             component_id < static_cast<int>(vec_var_n_components[variable_id]);
              ++component_id)
         {
             auto const global_component_id =
@@ -344,8 +351,8 @@ LocalToGlobalIndexMap::RowColumnIndices LocalToGlobalIndexMap::operator()(
                             _columns(mesh_item_id, component_id));
 }
 
-std::size_t
-LocalToGlobalIndexMap::getNumberOfElementDOF(std::size_t const mesh_item_id) const
+std::size_t LocalToGlobalIndexMap::getNumberOfElementDOF(
+    std::size_t const mesh_item_id) const
 {
     std::size_t ndof = 0;
 
@@ -357,8 +364,8 @@ LocalToGlobalIndexMap::getNumberOfElementDOF(std::size_t const mesh_item_id) con
     return ndof;
 }
 
-std::size_t
-LocalToGlobalIndexMap::getNumberOfElementComponents(std::size_t const mesh_item_id) const
+std::size_t LocalToGlobalIndexMap::getNumberOfElementComponents(
+    std::size_t const mesh_item_id) const
 {
     std::size_t n = 0;
     for (Table::Index c = 0; c < _rows.cols(); ++c)
@@ -377,7 +384,7 @@ std::vector<int> LocalToGlobalIndexMap::getElementVariableIDs(
     std::vector<int> vec;
     for (int i = 0; i < getNumberOfVariables(); i++)
     {
-        for (int j=0; j<getNumberOfVariableComponents(i); j++)
+        for (int j = 0; j < getNumberOfVariableComponents(i); j++)
         {
             auto comp_id = getGlobalComponent(i, j);
             if (!_rows(mesh_item_id, comp_id).empty())
@@ -450,8 +457,8 @@ std::ostream& operator<<(std::ostream& os, LocalToGlobalIndexMap const& map)
     std::size_t lines_printed = 0;
 
     os << "Rows of the local to global index map; " << map._rows.size()
-        << " rows\n";
-    for (std::size_t e=0; e<map.size(); ++e)
+       << " rows\n";
+    for (std::size_t e = 0; e < map.size(); ++e)
     {
         os << "== e " << e << " ==\n";
         for (int c = 0; c < map.getNumberOfGlobalComponents(); ++c)
@@ -460,7 +467,7 @@ std::ostream& operator<<(std::ostream& os, LocalToGlobalIndexMap const& map)
 
             os << "c" << c << " { ";
             std::copy(line.cbegin(), line.cend(),
-                std::ostream_iterator<std::size_t>(os, " "));
+                      std::ostream_iterator<std::size_t>(os, " "));
             os << " }\n";
         }
 
@@ -476,4 +483,4 @@ std::ostream& operator<<(std::ostream& os, LocalToGlobalIndexMap const& map)
 }
 #endif  // NDEBUG
 
-}   // namespace NumLib
+}  // namespace NumLib
diff --git a/NumLib/DOF/SimpleMatrixVectorProvider.cpp b/NumLib/DOF/SimpleMatrixVectorProvider.cpp
index 67fe93e71545deeaa79bc6d978f9f34e984d8e5b..9a0590099796995398ef789e296396c98f91def7 100644
--- a/NumLib/DOF/SimpleMatrixVectorProvider.cpp
+++ b/NumLib/DOF/SimpleMatrixVectorProvider.cpp
@@ -11,9 +11,9 @@
 #include "SimpleMatrixVectorProvider.h"
 
 #include <cassert>
-#include "BaseLib/Logging.h"
 
 #include "BaseLib/Error.h"
+#include "BaseLib/Logging.h"
 #include "MathLib/LinAlg/LinAlg.h"
 #include "MathLib/LinAlg/MatrixVectorTraits.h"
 
@@ -21,12 +21,10 @@ namespace LinAlg = MathLib::LinAlg;
 
 namespace detail
 {
-
-template<typename MatVec>
-MatVec*
-transfer(std::map<std::size_t, MatVec*>& from_unused,
-         std::map<MatVec*, std::size_t>& to_used,
-         typename std::map<std::size_t, MatVec*>::iterator it)
+template <typename MatVec>
+MatVec* transfer(std::map<std::size_t, MatVec*>& from_unused,
+                 std::map<MatVec*, std::size_t>& to_used,
+                 typename std::map<std::size_t, MatVec*>::iterator it)
 {
     auto const id = it->first;
     auto& ptr = it->second;
@@ -37,18 +35,17 @@ transfer(std::map<std::size_t, MatVec*>& from_unused,
     return res.first->first;
 }
 
-template<typename MatVec>
-void
-transfer(std::map<MatVec*, std::size_t>& from_used,
-         std::map<std::size_t, MatVec*>& to_unused,
-         typename std::map<MatVec*, std::size_t>::iterator it)
+template <typename MatVec>
+void transfer(std::map<MatVec*, std::size_t>& from_used,
+              std::map<std::size_t, MatVec*>& to_unused,
+              typename std::map<MatVec*, std::size_t>::iterator it)
 {
     auto& ptr = it->first;
     auto const id = it->second;
 
     auto res = to_unused.emplace(id, ptr);
     assert(res.second && "Emplacement failed.");
-    (void) res; // res unused if NDEBUG
+    (void)res;  // res unused if NDEBUG
     from_used.erase(it);
 }
 
@@ -56,17 +53,17 @@ transfer(std::map<MatVec*, std::size_t>& from_used,
 
 namespace NumLib
 {
-
-template<bool do_search, typename MatVec, typename... Args>
-std::pair<MatVec*, bool>
-SimpleMatrixVectorProvider::
-get_(std::size_t& id,
-     std::map<std::size_t, MatVec*>& unused_map,
-     std::map<MatVec*, std::size_t>& used_map,
-     Args&&... args)
+template <bool do_search, typename MatVec, typename... Args>
+std::pair<MatVec*, bool> SimpleMatrixVectorProvider::get_(
+    std::size_t& id,
+    std::map<std::size_t, MatVec*>& unused_map,
+    std::map<MatVec*, std::size_t>& used_map,
+    Args&&... args)
 {
-    if (id >= _next_id) {
-        OGS_FATAL("An obviously uninitialized id argument has been passed."
+    if (id >= _next_id)
+    {
+        OGS_FATAL(
+            "An obviously uninitialized id argument has been passed."
             " This might not be a serious error for the current implementation,"
             " but it might become one in the future."
             " Hence, I will abort now.");
@@ -83,57 +80,50 @@ get_(std::size_t& id,
 
     // not searched or not found, so create a new one
     id = _next_id++;
-    auto res = used_map.emplace(
-        MathLib::MatrixVectorTraits<MatVec>::newInstance(std::forward<Args>(args)...).release(),
-        id);
+    auto res =
+        used_map.emplace(MathLib::MatrixVectorTraits<MatVec>::newInstance(
+                             std::forward<Args>(args)...)
+                             .release(),
+                         id);
     assert(res.second && "Emplacement failed.");
-    return { res.first->first, true };
+    return {res.first->first, true};
 }
 
-template<bool do_search, typename... Args>
-std::pair<GlobalMatrix*, bool>
-SimpleMatrixVectorProvider::
-getMatrix_(std::size_t& id, Args&&... args)
+template <bool do_search, typename... Args>
+std::pair<GlobalMatrix*, bool> SimpleMatrixVectorProvider::getMatrix_(
+    std::size_t& id, Args&&... args)
 {
-    return get_<do_search>(id, _unused_matrices, _used_matrices, std::forward<Args>(args)...);
+    return get_<do_search>(
+        id, _unused_matrices, _used_matrices, std::forward<Args>(args)...);
 }
 
-
-GlobalMatrix&
-SimpleMatrixVectorProvider::
-getMatrix()
+GlobalMatrix& SimpleMatrixVectorProvider::getMatrix()
 {
     std::size_t id = 0u;
     return *getMatrix_<false>(id).first;
 }
 
-GlobalMatrix&
-SimpleMatrixVectorProvider::
-getMatrix(std::size_t& id)
+GlobalMatrix& SimpleMatrixVectorProvider::getMatrix(std::size_t& id)
 {
     return *getMatrix_<true>(id).first;
 }
 
-GlobalMatrix&
-SimpleMatrixVectorProvider::
-getMatrix(MathLib::MatrixSpecifications const& ms)
+GlobalMatrix& SimpleMatrixVectorProvider::getMatrix(
+    MathLib::MatrixSpecifications const& ms)
 {
     std::size_t id = 0u;
     return *getMatrix_<false>(id, ms).first;
     // TODO assert that the returned object always is of the right size
 }
 
-GlobalMatrix&
-SimpleMatrixVectorProvider::
-getMatrix(MathLib::MatrixSpecifications const& ms, std::size_t& id)
+GlobalMatrix& SimpleMatrixVectorProvider::getMatrix(
+    MathLib::MatrixSpecifications const& ms, std::size_t& id)
 {
     return *getMatrix_<true>(id, ms).first;
     // TODO assert that the returned object always is of the right size
 }
 
-GlobalMatrix&
-SimpleMatrixVectorProvider::
-getMatrix(GlobalMatrix const& A)
+GlobalMatrix& SimpleMatrixVectorProvider::getMatrix(GlobalMatrix const& A)
 {
     std::size_t id = 0u;
     auto const& res = getMatrix_<false>(id, A);
@@ -144,9 +134,8 @@ getMatrix(GlobalMatrix const& A)
     return *res.first;
 }
 
-GlobalMatrix&
-SimpleMatrixVectorProvider::
-getMatrix(GlobalMatrix const& A, std::size_t& id)
+GlobalMatrix& SimpleMatrixVectorProvider::getMatrix(GlobalMatrix const& A,
+                                                    std::size_t& id)
 {
     auto const& res = getMatrix_<true>(id, A);
     if (!res.second)
@@ -156,53 +145,50 @@ getMatrix(GlobalMatrix const& A, std::size_t& id)
     return *res.first;
 }
 
-void
-SimpleMatrixVectorProvider::
-releaseMatrix(GlobalMatrix const& A)
+void SimpleMatrixVectorProvider::releaseMatrix(GlobalMatrix const& A)
 {
     auto it = _used_matrices.find(const_cast<GlobalMatrix*>(&A));
-    if (it == _used_matrices.end()) {
-        OGS_FATAL("The given matrix has not been found. Cannot release it. Aborting.");
-    } else {
+    if (it == _used_matrices.end())
+    {
+        OGS_FATAL(
+            "The given matrix has not been found. Cannot release it. "
+            "Aborting.");
+    }
+    else
+    {
         ::detail::transfer(_used_matrices, _unused_matrices, it);
     }
 }
 
-template<bool do_search, typename... Args>
-std::pair<GlobalVector*, bool>
-SimpleMatrixVectorProvider::
-getVector_(std::size_t& id, Args&&... args)
+template <bool do_search, typename... Args>
+std::pair<GlobalVector*, bool> SimpleMatrixVectorProvider::getVector_(
+    std::size_t& id, Args&&... args)
 {
-    return get_<do_search>(id, _unused_vectors, _used_vectors, std::forward<Args>(args)...);
+    return get_<do_search>(
+        id, _unused_vectors, _used_vectors, std::forward<Args>(args)...);
 }
 
-GlobalVector&
-SimpleMatrixVectorProvider::
-getVector(std::size_t& id)
+GlobalVector& SimpleMatrixVectorProvider::getVector(std::size_t& id)
 {
     return *getVector_<true>(id).first;
 }
 
-GlobalVector&
-SimpleMatrixVectorProvider::
-getVector(MathLib::MatrixSpecifications const& ms)
+GlobalVector& SimpleMatrixVectorProvider::getVector(
+    MathLib::MatrixSpecifications const& ms)
 {
     std::size_t id = 0u;
     return *getVector_<false>(id, ms).first;
     // TODO assert that the returned object always is of the right size
 }
 
-GlobalVector&
-SimpleMatrixVectorProvider::
-getVector(MathLib::MatrixSpecifications const& ms, std::size_t& id)
+GlobalVector& SimpleMatrixVectorProvider::getVector(
+    MathLib::MatrixSpecifications const& ms, std::size_t& id)
 {
     return *getVector_<true>(id, ms).first;
     // TODO assert that the returned object always is of the right size
 }
 
-GlobalVector&
-SimpleMatrixVectorProvider::
-getVector(GlobalVector const& x)
+GlobalVector& SimpleMatrixVectorProvider::getVector(GlobalVector const& x)
 {
     std::size_t id = 0u;
     auto const& res = getVector_<false>(id, x);
@@ -213,9 +199,8 @@ getVector(GlobalVector const& x)
     return *res.first;
 }
 
-GlobalVector&
-SimpleMatrixVectorProvider::
-getVector(GlobalVector const& x, std::size_t& id)
+GlobalVector& SimpleMatrixVectorProvider::getVector(GlobalVector const& x,
+                                                    std::size_t& id)
 {
     auto const& res = getVector_<true>(id, x);
     if (!res.second)
@@ -225,20 +210,22 @@ getVector(GlobalVector const& x, std::size_t& id)
     return *res.first;
 }
 
-void
-SimpleMatrixVectorProvider::
-releaseVector(GlobalVector const& x)
+void SimpleMatrixVectorProvider::releaseVector(GlobalVector const& x)
 {
     auto it = _used_vectors.find(const_cast<GlobalVector*>(&x));
-    if (it == _used_vectors.end()) {
-        OGS_FATAL("The given vector has not been found. Cannot release it. Aborting.");
-    } else {
+    if (it == _used_vectors.end())
+    {
+        OGS_FATAL(
+            "The given vector has not been found. Cannot release it. "
+            "Aborting.");
+    }
+    else
+    {
         ::detail::transfer(_used_vectors, _unused_vectors, it);
     }
 }
 
-SimpleMatrixVectorProvider::
-~SimpleMatrixVectorProvider()
+SimpleMatrixVectorProvider::~SimpleMatrixVectorProvider()
 {
     if (!_used_matrices.empty())
     {
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
index bf63356160281d163516b2272640d6f01c43d1fd..a742023871cb766c07834961b6d6e9d82979df8d 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
@@ -11,14 +11,14 @@
 #include "LocalLinearLeastSquaresExtrapolator.h"
 
 #include <Eigen/SVD>
-#include "BaseLib/Logging.h"
 
+#include "BaseLib/Logging.h"
+#include "ExtrapolatableElementCollection.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "MathLib/LinAlg/LinAlg.h"
 #include "MathLib/LinAlg/MatrixVectorTraits.h"
 #include "NumLib/Assembler/SerialExecutor.h"
 #include "NumLib/Function/Interpolation.h"
-#include "ExtrapolatableElementCollection.h"
 
 namespace NumLib
 {
diff --git a/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp b/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp
index 8d79c6cac55ee9a9b4353ab192ec8c83fcee39b2..b165dad5af2aea9a4bedc9bc11cea362258ace7c 100644
--- a/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp
+++ b/NumLib/Fem/CoordinatesMapping/NaturalCoordinatesMapping.cpp
@@ -16,9 +16,7 @@
 #endif  // NDEBUG
 
 #include "BaseLib/Error.h"
-
 #include "MeshLib/ElementCoordinatesMappingLocal.h"
-#include "MeshLib/Elements/TemplateElement.h"
 #include "MeshLib/Elements/HexRule20.h"
 #include "MeshLib/Elements/HexRule8.h"
 #include "MeshLib/Elements/LineRule2.h"
@@ -31,38 +29,38 @@
 #include "MeshLib/Elements/QuadRule4.h"
 #include "MeshLib/Elements/QuadRule8.h"
 #include "MeshLib/Elements/QuadRule9.h"
+#include "MeshLib/Elements/TemplateElement.h"
 #include "MeshLib/Elements/TetRule10.h"
 #include "MeshLib/Elements/TetRule4.h"
 #include "MeshLib/Elements/TriRule3.h"
 #include "MeshLib/Elements/TriRule6.h"
-
-#include "NumLib/Fem/ShapeFunction/ShapePoint1.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
+#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
 #include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
 #include "NumLib/Fem/ShapeFunction/ShapeLine3.h"
-#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
-#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
+#include "NumLib/Fem/ShapeFunction/ShapePoint1.h"
+#include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
+#include "NumLib/Fem/ShapeFunction/ShapePrism6.h"
+#include "NumLib/Fem/ShapeFunction/ShapePyra13.h"
+#include "NumLib/Fem/ShapeFunction/ShapePyra5.h"
 #include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
 #include "NumLib/Fem/ShapeFunction/ShapeQuad8.h"
 #include "NumLib/Fem/ShapeFunction/ShapeQuad9.h"
-#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
-#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
-#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
 #include "NumLib/Fem/ShapeFunction/ShapeTet10.h"
-#include "NumLib/Fem/ShapeFunction/ShapePrism6.h"
-#include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
-#include "NumLib/Fem/ShapeFunction/ShapePyra5.h"
-#include "NumLib/Fem/ShapeFunction/ShapePyra13.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
+#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
-
 #include "ShapeMatrices.h"
 
 namespace NumLib
 {
-
 namespace detail
 {
-
-template <ShapeMatrixType FIELD_TYPE> struct FieldType {};
+template <ShapeMatrixType FIELD_TYPE>
+struct FieldType
+{
+};
 
 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
 inline void computeMappingMatrices(
@@ -107,8 +105,9 @@ static void checkJacobianDeterminant(const double detJ,
         std::cerr << element << "\n";
 #endif  // NDEBUG
         OGS_FATAL(
-                "Please check whether the node numbering of the element is correct,"
-                "or additional elements (like boundary elements) are still present in the mesh.");
+            "Please check whether the node numbering of the element is correct,"
+            "or additional elements (like boundary elements) are still present "
+            "in the mesh.");
     }
 
     if (detJ == 0)
@@ -180,10 +179,18 @@ inline void computeMappingMatrices(
     T_SHAPE_MATRICES& shapemat,
     FieldType<ShapeMatrixType::N_J> /*unused*/)
 {
-    computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
-        (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::N>());
-    computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
-        (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
+    computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>(
+        ele,
+        natural_pt,
+        ele_local_coord,
+        shapemat,
+        FieldType<ShapeMatrixType::N>());
+    computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>(
+        ele,
+        natural_pt,
+        ele_local_coord,
+        shapemat,
+        FieldType<ShapeMatrixType::DNDR_J>());
 }
 
 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
@@ -236,10 +243,18 @@ inline void computeMappingMatrices(
     T_SHAPE_MATRICES& shapemat,
     FieldType<ShapeMatrixType::ALL> /*unused*/)
 {
-    computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
-        (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::N>());
-    computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
-        (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDX>());
+    computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>(
+        ele,
+        natural_pt,
+        ele_local_coord,
+        shapemat,
+        FieldType<ShapeMatrixType::N>());
+    computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>(
+        ele,
+        natural_pt,
+        ele_local_coord,
+        shapemat,
+        FieldType<ShapeMatrixType::DNDX>());
 }
 
 template <class T_MESH_ELEMENT,
@@ -251,17 +266,16 @@ void naturalCoordinatesMappingComputeShapeMatrices(const T_MESH_ELEMENT& ele,
                                                    T_SHAPE_MATRICES& shapemat,
                                                    const unsigned global_dim)
 {
-    const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(ele, global_dim);
+    const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(ele,
+                                                                  global_dim);
 
-    detail::computeMappingMatrices<
-        T_MESH_ELEMENT,
-        T_SHAPE_FUNC,
-        T_SHAPE_MATRICES>
-            (ele,
-             natural_pt,
-             ele_local_coord,
-             shapemat,
-             detail::FieldType<T_SHAPE_MATRIX_TYPE>());
+    detail::
+        computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>(
+            ele,
+            natural_pt,
+            ele_local_coord,
+            shapemat,
+            detail::FieldType<T_SHAPE_MATRIX_TYPE>());
 }
 
 #define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART(        \
diff --git a/NumLib/Fem/ShapeFunction/ShapeStaticConsts.cpp b/NumLib/Fem/ShapeFunction/ShapeStaticConsts.cpp
index 3956551c193d942c0f98bee0c1264389b6358dc7..a25830b9ab499bb6949c44c8ce2da3941ce82a8e 100644
--- a/NumLib/Fem/ShapeFunction/ShapeStaticConsts.cpp
+++ b/NumLib/Fem/ShapeFunction/ShapeStaticConsts.cpp
@@ -47,7 +47,6 @@
 
 namespace NumLib
 {
-
 const unsigned ShapeHex20::DIM;
 const unsigned ShapeHex20::NPOINTS;
 const unsigned ShapeHex8::DIM;
diff --git a/NumLib/ODESolver/ConvergenceCriterion.cpp b/NumLib/ODESolver/ConvergenceCriterion.cpp
index a09c416f98d6923b8a02cc1f8379fd950bee9eaa..ab09bec442e21f301c25306469be4d74abffed7a 100644
--- a/NumLib/ODESolver/ConvergenceCriterion.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterion.cpp
@@ -9,13 +9,13 @@
  */
 
 #include "ConvergenceCriterion.h"
+
 #include "BaseLib/ConfigTree.h"
 #include "BaseLib/Error.h"
-
 #include "ConvergenceCriterionDeltaX.h"
-#include "ConvergenceCriterionResidual.h"
 #include "ConvergenceCriterionPerComponentDeltaX.h"
 #include "ConvergenceCriterionPerComponentResidual.h"
+#include "ConvergenceCriterionResidual.h"
 
 namespace NumLib
 {
@@ -25,7 +25,8 @@ std::unique_ptr<ConvergenceCriterion> createConvergenceCriterion(
     //! \ogs_file_param{prj__time_loop__processes__process__convergence_criterion__type}
     auto const type = config.peekConfigParameter<std::string>("type");
 
-    if (type == "DeltaX") {
+    if (type == "DeltaX")
+    {
         return createConvergenceCriterionDeltaX(config);
     }
     if (type == "Residual")
diff --git a/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp b/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp
index 24074d10e6ac3465178cf615c9f5e4f517d38821..20f2e03b938f19cfbd0d5e7c1c0a41084b084828 100644
--- a/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp
@@ -9,9 +9,9 @@
  */
 
 #include "ConvergenceCriterionDeltaX.h"
-#include "BaseLib/Logging.h"
 
 #include "BaseLib/ConfigTree.h"
+#include "BaseLib/Logging.h"
 #include "MathLib/LinAlg/LinAlg.h"
 
 namespace NumLib
@@ -46,10 +46,12 @@ void ConvergenceCriterionDeltaX::checkDeltaX(const GlobalVector& minus_delta_x,
     bool satisfied_abs = false;
     bool satisfied_rel = false;
 
-    if (_abstol) {
+    if (_abstol)
+    {
         satisfied_abs = error_dx < *_abstol;
     }
-    if (_reltol) {
+    if (_reltol)
+    {
         satisfied_rel = checkRelativeTolerance(*_reltol, error_dx, norm_x);
     }
 
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp
index 56b6724eab3c6f9d5c0348f2e6c64e10effc58b6..25f875333dbfe984d825bf6ad140c1f81ff373e5 100644
--- a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp
@@ -9,12 +9,13 @@
  */
 
 #include "ConvergenceCriterionPerComponentDeltaX.h"
+
 #include <limits>
-#include "BaseLib/Logging.h"
 
 #include "BaseLib/ConfigTree.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "BaseLib/Logging.h"
 #include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
 
 namespace NumLib
 {
@@ -112,9 +113,12 @@ createConvergenceCriterionPerComponentDeltaX(const BaseLib::ConfigTree& config)
             "At least one of absolute or relative tolerance has to be "
             "specified.");
     }
-    if (!abstols) {
+    if (!abstols)
+    {
         abstols = std::vector<double>(reltols->size());
-    } else if (!reltols) {
+    }
+    else if (!reltols)
+    {
         reltols = std::vector<double>(abstols->size());
     }
 
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
index db24a3b2c57b41ffd95be90aa7ea85fd75fcc37e..c389b5a585fc02200ce5aeaf50fd766d6341db74 100644
--- a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
@@ -9,9 +9,9 @@
  */
 
 #include "ConvergenceCriterionPerComponentResidual.h"
-#include "BaseLib/Logging.h"
 
 #include "BaseLib/ConfigTree.h"
+#include "BaseLib/Logging.h"
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
 
@@ -40,7 +40,6 @@ ConvergenceCriterionPerComponentResidual::
     }
 }
 
-
 void ConvergenceCriterionPerComponentResidual::checkDeltaX(
     const GlobalVector& minus_delta_x, GlobalVector const& x)
 {
@@ -67,7 +66,6 @@ void ConvergenceCriterionPerComponentResidual::checkDeltaX(
     }
 }
 
-
 void ConvergenceCriterionPerComponentResidual::checkResidual(
     const GlobalVector& residual)
 {
@@ -85,14 +83,17 @@ void ConvergenceCriterionPerComponentResidual::checkResidual(
          ++global_component)
     {
         // TODO short cut if tol <= 0.0
-        auto norm_res = norm(residual, global_component, _norm_type,
-                             *_dof_table, *_mesh);
+        auto norm_res =
+            norm(residual, global_component, _norm_type, *_dof_table, *_mesh);
 
-        if (_is_first_iteration) {
+        if (_is_first_iteration)
+        {
             INFO("Convergence criterion, component {:d}: |r0|={:.4e}",
                  global_component, norm_res);
             _residual_norms_0[global_component] = norm_res;
-        } else {
+        }
+        else
+        {
             auto const norm_res0 = _residual_norms_0[global_component];
             INFO(
                 "Convergence criterion, component {:d}: |r|={:.4e}, "
@@ -151,9 +152,12 @@ createConvergenceCriterionPerComponentResidual(
             "At least one of absolute or relative tolerance has to be "
             "specified.");
     }
-    if (!abstols) {
+    if (!abstols)
+    {
         abstols = std::vector<double>(reltols->size());
-    } else if (!reltols) {
+    }
+    else if (!reltols)
+    {
         reltols = std::vector<double>(abstols->size());
     }
 
diff --git a/NumLib/ODESolver/ConvergenceCriterionResidual.cpp b/NumLib/ODESolver/ConvergenceCriterionResidual.cpp
index 08de36c59ab5dd72bb8b599ec816627656a9846f..125d62eceaf0f463c35c9ed03e8b5656562e6a67 100644
--- a/NumLib/ODESolver/ConvergenceCriterionResidual.cpp
+++ b/NumLib/ODESolver/ConvergenceCriterionResidual.cpp
@@ -9,9 +9,9 @@
  */
 
 #include "ConvergenceCriterionResidual.h"
-#include "BaseLib/Logging.h"
 
 #include "BaseLib/ConfigTree.h"
+#include "BaseLib/Logging.h"
 #include "MathLib/LinAlg/LinAlg.h"
 
 namespace NumLib
diff --git a/NumLib/ODESolver/MatrixTranslator.cpp b/NumLib/ODESolver/MatrixTranslator.cpp
index 76d1cebe2bd3e2a309dd61df57c832ab5156c1a8..81f68482cbe8c15ece7390bd5368bcfbc9b9686c 100644
--- a/NumLib/ODESolver/MatrixTranslator.cpp
+++ b/NumLib/ODESolver/MatrixTranslator.cpp
@@ -52,7 +52,7 @@ void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
 
     // res = M * x_dot + K * x_curr - b
     LinAlg::matMult(M, xdot, res);  // the local vector x_dot seems to be
-                                  // necessary because of this multiplication
+                                    // necessary because of this multiplication
     LinAlg::matMultAdd(K, x_curr, res, res);
     LinAlg::axpy(res, -1.0, b);
 }
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index 6a8ffcd420afcdd6813fc01fbe47eb3e62127571..760901b22f43d158aeb38a793e1a8f7d378a2e3a 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -56,10 +56,8 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
     namespace LinAlg = MathLib::LinAlg;
     auto& sys = *_equation_system;
 
-    auto& A =
-        NumLib::GlobalMatrixProvider::provider.getMatrix(_A_id);
-    auto& rhs = NumLib::GlobalVectorProvider::provider.getVector(
-        _rhs_id);
+    auto& A = NumLib::GlobalMatrixProvider::provider.getMatrix(_A_id);
+    auto& rhs = NumLib::GlobalVectorProvider::provider.getVector(_rhs_id);
 
     std::vector<GlobalVector*> x_new{x};
     x_new[process_id] =
@@ -71,8 +69,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
     _convergence_criterion->preFirstIteration();
 
     int iteration = 1;
-    for (; iteration <= _maxiter;
-         ++iteration, _convergence_criterion->reset())
+    for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
     {
         BaseLib::RunTime timer_dirichlet;
         double time_dirichlet = 0.0;
@@ -105,10 +102,11 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
         time_dirichlet += timer_dirichlet.elapsed();
         INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
 
-        if (!sys.isLinear() && _convergence_criterion->hasResidualCheck()) {
+        if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
+        {
             GlobalVector res;
             LinAlg::matMult(A, *x_new[process_id], res);  // res = A * x_new
-            LinAlg::axpy(res, -1.0, rhs);   // res -= rhs
+            LinAlg::axpy(res, -1.0, rhs);                 // res -= rhs
             _convergence_criterion->checkResidual(res);
         }
 
@@ -162,10 +160,14 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
             break;
         }
 
-        if (sys.isLinear()) {
+        if (sys.isLinear())
+        {
             error_norms_met = true;
-        } else {
-            if (_convergence_criterion->hasDeltaXCheck()) {
+        }
+        else
+        {
+            if (_convergence_criterion->hasDeltaXCheck())
+            {
                 GlobalVector minus_delta_x(*x[process_id]);
                 LinAlg::axpy(minus_delta_x, -1.0,
                              *x_new[process_id]);  // minus_delta_x = x - x_new
@@ -234,13 +236,10 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
     namespace LinAlg = MathLib::LinAlg;
     auto& sys = *_equation_system;
 
-    auto& res = NumLib::GlobalVectorProvider::provider.getVector(
-        _res_id);
+    auto& res = NumLib::GlobalVectorProvider::provider.getVector(_res_id);
     auto& minus_delta_x =
-        NumLib::GlobalVectorProvider::provider.getVector(
-            _minus_delta_x_id);
-    auto& J =
-        NumLib::GlobalMatrixProvider::provider.getMatrix(_J_id);
+        NumLib::GlobalVectorProvider::provider.getVector(_minus_delta_x_id);
+    auto& J = NumLib::GlobalMatrixProvider::provider.getMatrix(_J_id);
 
     bool error_norms_met = false;
 
@@ -251,8 +250,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
     _convergence_criterion->preFirstIteration();
 
     int iteration = 1;
-    for (; iteration <= _maxiter;
-         ++iteration, _convergence_criterion->reset())
+    for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
     {
         BaseLib::RunTime timer_dirichlet;
         double time_dirichlet = 0.0;
@@ -361,10 +359,14 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
             break;
         }
 
-        if (sys.isLinear()) {
+        if (sys.isLinear())
+        {
             error_norms_met = true;
-        } else {
-            if (_convergence_criterion->hasDeltaXCheck()) {
+        }
+        else
+        {
+            if (_convergence_criterion->hasDeltaXCheck())
+            {
                 // Note: x contains the new solution!
                 _convergence_criterion->checkDeltaX(minus_delta_x,
                                                     *x[process_id]);
@@ -398,8 +400,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
 
     NumLib::GlobalMatrixProvider::provider.releaseMatrix(J);
     NumLib::GlobalVectorProvider::provider.releaseVector(res);
-    NumLib::GlobalVectorProvider::provider.releaseVector(
-        minus_delta_x);
+    NumLib::GlobalVectorProvider::provider.releaseVector(minus_delta_x);
 
     return {error_norms_met, iteration};
 }
@@ -413,7 +414,8 @@ createNonlinearSolver(GlobalLinearSolver& linear_solver,
     //! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__max_iter}
     auto const max_iter = config.getConfigParameter<int>("max_iter");
 
-    if (type == "Picard") {
+    if (type == "Picard")
+    {
         auto const tag = NonlinearSolverTag::Picard;
         using ConcreteNLS = NonlinearSolver<tag>;
         return std::make_pair(
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index d28ea348a1c45f9a2e17e8b815ddb70edc619b3d..3d941cd584d6ada88e5e5b9fc4c51d67439af991 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -12,8 +12,8 @@
 
 #include "MathLib/LinAlg/ApplyKnownSolution.h"
 #include "MathLib/LinAlg/UnifiedMatrixSetters.h"
-#include "NumLib/IndexValueVector.h"
 #include "NumLib/Exceptions.h"
+#include "NumLib/IndexValueVector.h"
 
 namespace detail
 {
@@ -87,7 +87,8 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     _xdot_ids.resize(x_new_timestep.size());
     for (std::size_t i = 0; i < xdot.size(); i++)
     {
-        xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector(_xdot_ids[i]);
+        xdot[i] =
+            &NumLib::GlobalVectorProvider::provider.getVector(_xdot_ids[i]);
         _time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
     }
 
@@ -225,7 +226,8 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
 
     for (std::size_t i = 0; i < xdot.size(); i++)
     {
-        xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector(_xdot_ids[i]);
+        xdot[i] =
+            &NumLib::GlobalVectorProvider::provider.getVector(_xdot_ids[i]);
         _time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
     }
 
diff --git a/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp b/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp
index ed5a03ccbefe1624ee679ad365d719f398fc9651..c079d88b56479bdaba02102e678f2e8830603316 100644
--- a/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp
+++ b/NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.cpp
@@ -13,7 +13,6 @@
 
 #include "BaseLib/Algorithm.h"
 #include "BaseLib/ConfigTree.h"
-
 #include "EvolutionaryPIDcontroller.h"
 #include "TimeStepAlgorithm.h"
 
diff --git a/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp
index d4ea86361e44624f4a87f4f5271365afbc87960f..32c1d34142614071fd35e99611605524980cb9af 100644
--- a/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.cpp
@@ -10,11 +10,11 @@
  */
 
 #include "CreateFixedTimeStepping.h"
+
 #include <string>
 
 #include "BaseLib/ConfigTree.h"
 #include "BaseLib/Error.h"
-
 #include "FixedTimeStepping.h"
 #include "TimeStepAlgorithm.h"
 
diff --git a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
index 89463aef6d662aa140c713853e0fb730bda064d4..05feeda9b6acc8ce2260354e91c3f78748b58b89 100644
--- a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
+++ b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
@@ -14,9 +14,9 @@
 #include <functional>
 #include <limits>
 #include <vector>
-#include "BaseLib/Logging.h"
 
 #include "BaseLib/Algorithm.h"
+#include "BaseLib/Logging.h"
 
 namespace NumLib
 {
diff --git a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
index e1a1a463c976d5750cc60d4cbf1ee1c9d9fb8b68..038c84e1249a1b7508ec4e9129deacbaa7be75ae 100644
--- a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
@@ -13,9 +13,9 @@
 #include "FixedTimeStepping.h"
 
 #include <algorithm>
-#include <numeric>
-#include <limits>
 #include <cassert>
+#include <limits>
+#include <numeric>
 
 namespace NumLib
 {
diff --git a/NumLib/TimeStepping/CreateTimeStepper.cpp b/NumLib/TimeStepping/CreateTimeStepper.cpp
index 384cda3c5a77eb8b09bc97e078924c1a02617d28..4074160c5695e2ba17824db00635aa05642c31ba 100644
--- a/NumLib/TimeStepping/CreateTimeStepper.cpp
+++ b/NumLib/TimeStepping/CreateTimeStepper.cpp
@@ -16,7 +16,6 @@
 
 #include "BaseLib/ConfigTree.h"
 #include "BaseLib/Error.h"
-
 #include "NumLib/TimeStepping/Algorithms/CreateEvolutionaryPIDcontroller.h"
 #include "NumLib/TimeStepping/Algorithms/CreateFixedTimeStepping.h"
 #include "NumLib/TimeStepping/Algorithms/CreateIterationNumberBasedTimeStepping.h"