diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/Utils/CreateLocalAssemblers.h b/ProcessLib/BoundaryConditionAndSourceTerm/Utils/CreateLocalAssemblers.h
index d6299bbfea2bc4237dbe67bee3346eb3da9750e4..37e027f7c5f2a7bb994356cafc831bc62ede82ed 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/Utils/CreateLocalAssemblers.h
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/Utils/CreateLocalAssemblers.h
@@ -12,7 +12,7 @@
 #include <vector>
 
 #include "BaseLib/Logging.h"
-#include "LocalDataInitializer.h"
+#include "LocalAssemblerFactory.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
 
 namespace ProcessLib
@@ -34,21 +34,20 @@ void createLocalAssemblers(
     static_assert(
         GlobalDim == 1 || GlobalDim == 2 || GlobalDim == 3,
         "Meshes with dimension greater than three are not supported.");
-    // Shape matrices initializer
-    using LocalDataInitializer =
-        LocalDataInitializer<LocalAssemblerInterface,
-                             LocalAssemblerImplementation, GlobalDim,
-                             ExtraCtorArgs...>;
+
+    using LocalAssemblerFactory =
+        LocalAssemblerFactory<LocalAssemblerInterface,
+                              LocalAssemblerImplementation, GlobalDim,
+                              ExtraCtorArgs...>;
 
     DBUG("Create local assemblers.");
-    // Populate the vector of local assemblers.
-    local_assemblers.resize(mesh_elements.size());
 
-    LocalDataInitializer initializer(dof_table, shapefunction_order);
+    LocalAssemblerFactory factory(dof_table, shapefunction_order);
+    local_assemblers.resize(mesh_elements.size());
 
     DBUG("Calling local assembler builder for all mesh elements.");
     GlobalExecutor::transformDereferenced(
-        initializer, mesh_elements, local_assemblers,
+        factory, mesh_elements, local_assemblers,
         std::forward<ExtraCtorArgs>(extra_ctor_args)...);
 }
 }  // namespace detail
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/Utils/LocalAssemblerFactory.h b/ProcessLib/BoundaryConditionAndSourceTerm/Utils/LocalAssemblerFactory.h
new file mode 100644
index 0000000000000000000000000000000000000000..3f0fd58d6e71216fedc35cb4889745a1e6ac11bb
--- /dev/null
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/Utils/LocalAssemblerFactory.h
@@ -0,0 +1,116 @@
+/**
+ * \file
+ * \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
+
+#include "ProcessLib/Utils/EnabledElements.h"
+#include "ProcessLib/Utils/GenericLocalAssemblerFactory.h"
+
+namespace ProcessLib
+{
+namespace BoundaryConditionAndSourceTerm
+{
+template <typename LocalAssemblerInterface,
+          template <typename, typename, int>
+          class LocalAssemblerImplementation,
+          int GlobalDim,
+          typename... ConstructorArgs>
+class LocalAssemblerFactory final
+    : public ProcessLib::GenericLocalAssemblerFactory<LocalAssemblerInterface,
+                                                      ConstructorArgs...>
+{
+    using Base =
+        ProcessLib::GenericLocalAssemblerFactory<LocalAssemblerInterface,
+                                                 ConstructorArgs...>;
+
+    template <typename ShapeFunction>
+    using LocAsmBuilderFactory =
+        ProcessLib::LocalAssemblerBuilderFactory<ShapeFunction,
+                                                 LocalAssemblerInterface,
+                                                 LocalAssemblerImplementation,
+                                                 GlobalDim,
+                                                 ConstructorArgs...>;
+
+    struct HasSuitableDimension
+    {
+        template <typename ElementTraits>
+        constexpr bool operator()(ElementTraits*) const
+        {
+            return GlobalDim >= ElementTraits::ShapeFunction::DIM;
+        }
+    };
+
+    struct Is2ndOrderElementOfSuitableDimension
+    {
+        template <typename ElementTraits>
+        constexpr bool operator()(ElementTraits*) const
+        {
+            if constexpr (GlobalDim < ElementTraits::ShapeFunction::DIM)
+            {
+                return false;
+            }
+            return ElementTraits::ShapeFunction::ORDER == 2 ||
+                   // points are needed for 2nd order, too
+                   std::is_same_v<MeshLib::Point,
+                                  typename ElementTraits::Element>;
+        }
+    };
+
+public:
+    LocalAssemblerFactory(NumLib::LocalToGlobalIndexMap const& dof_table,
+                          const unsigned shapefunction_order)
+        : Base(dof_table)
+    {
+        if (shapefunction_order < 1 || 2 < shapefunction_order)
+        {
+            OGS_FATAL("The given shape function order {:d} is not supported",
+                      shapefunction_order);
+        }
+
+        if (shapefunction_order == 1)
+        {
+            // 1st order is enabled on all elements with suitable dimension
+            using EnabledElementTraits =
+                decltype(BaseLib::TMP::filter<EnabledElementTraitsLagrange>(
+                    std::declval<HasSuitableDimension>()));
+
+            BaseLib::TMP::foreach<EnabledElementTraits>(
+                [this]<typename ET>(ET*)
+                {
+                    using MeshElement = typename ET::Element;
+                    // this will use linear shape functions on higher order
+                    // elements and the linear shape function on linear elements
+                    using LowerOrderShapeFunction =
+                        typename ET::LowerOrderShapeFunction;
+                    Base::_builders[std::type_index(typeid(MeshElement))] =
+                        LocAsmBuilderFactory<LowerOrderShapeFunction>::create();
+                });
+        }
+        else if (shapefunction_order == 2)
+        {
+            // 2nd order only on 2nd order elements
+            using EnabledElementTraits =
+                decltype(BaseLib::TMP::filter<EnabledElementTraitsLagrange>(
+                    std::declval<Is2ndOrderElementOfSuitableDimension>()));
+
+            BaseLib::TMP::foreach<EnabledElementTraits>(
+                [this]<typename ET>(ET*)
+                {
+                    using MeshElement = typename ET::Element;
+                    using ShapeFunction2ndOrder = typename ET::ShapeFunction;
+                    Base::_builders[std::type_index(typeid(MeshElement))] =
+                        LocAsmBuilderFactory<ShapeFunction2ndOrder>::create();
+                });
+        }
+    }
+};
+
+}  // namespace BoundaryConditionAndSourceTerm
+}  // namespace ProcessLib
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/Utils/LocalDataInitializer.h b/ProcessLib/BoundaryConditionAndSourceTerm/Utils/LocalDataInitializer.h
deleted file mode 100644
index 204b949d2ca33e8293303ad4e05c25b2c1043ad4..0000000000000000000000000000000000000000
--- a/ProcessLib/BoundaryConditionAndSourceTerm/Utils/LocalDataInitializer.h
+++ /dev/null
@@ -1,180 +0,0 @@
-/**
- * \file
- * \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
-
-#include <functional>
-#include <memory>
-#include <type_traits>
-#include <typeindex>
-#include <typeinfo>
-#include <unordered_map>
-
-#include "MeshLib/Elements/Elements.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
-#include "ProcessLib/Utils/EnabledElements.h"
-
-namespace ProcessLib
-{
-namespace BoundaryConditionAndSourceTerm
-{
-/// The LocalDataInitializer is a functor creating a local assembler data with
-/// corresponding to the mesh element type shape functions and calling
-/// initialization of the new local assembler data.
-/// For example for MeshLib::Quad a local assembler data with template argument
-/// NumLib::ShapeQuad4 is created.
-template <typename LocalAssemblerInterface,
-          template <typename, typename, int> class LocalAssemblerData,
-          int GlobalDim, typename... ConstructorArgs>
-class LocalDataInitializer final
-{
-    struct HasSuitableDimension
-    {
-        template <typename ElementTraits>
-        constexpr bool operator()(ElementTraits*) const
-        {
-            return GlobalDim >= ElementTraits::ShapeFunction::DIM;
-        }
-    };
-
-    struct Is2ndOrderElementOfSuitableDimension
-    {
-        template <typename ElementTraits>
-        constexpr bool operator()(ElementTraits*) const
-        {
-            if constexpr (GlobalDim < ElementTraits::ShapeFunction::DIM)
-            {
-                return false;
-            }
-            return ElementTraits::ShapeFunction::ORDER == 2 ||
-                   // points are needed for 2nd order, too
-                   std::is_same_v<MeshLib::Point,
-                                  typename ElementTraits::Element>;
-        }
-    };
-
-public:
-    using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
-
-    LocalDataInitializer(NumLib::LocalToGlobalIndexMap const& dof_table,
-                         const unsigned shapefunction_order)
-        : _dof_table(dof_table)
-    {
-        if (shapefunction_order < 1 || 2 < shapefunction_order)
-        {
-            OGS_FATAL("The given shape function order {:d} is not supported",
-                      shapefunction_order);
-        }
-
-        if (shapefunction_order == 1)
-        {
-            // 1st order is enabled on all elements with suitable dimension
-            using EnabledElementTraits =
-                decltype(BaseLib::TMP::filter<EnabledElementTraitsLagrange>(
-                    std::declval<HasSuitableDimension>()));
-
-            BaseLib::TMP::foreach<EnabledElementTraits>(
-                [this]<typename ET>(ET*)
-                {
-                    using MeshElement = typename ET::Element;
-                    // this will use linear shape functions on higher order
-                    // elements and the linear shape function on linear elements
-                    using LowerOrderShapeFunction =
-                        typename ET::LowerOrderShapeFunction;
-                    _builder[std::type_index(typeid(MeshElement))] =
-                        makeLocalAssemblerBuilder<LowerOrderShapeFunction>();
-                });
-        }
-        else if (shapefunction_order == 2)
-        {
-            // 2nd order only on 2nd order elements
-            using EnabledElementTraits =
-                decltype(BaseLib::TMP::filter<EnabledElementTraitsLagrange>(
-                    std::declval<Is2ndOrderElementOfSuitableDimension>()));
-
-            BaseLib::TMP::foreach<EnabledElementTraits>(
-                [this]<typename ET>(ET*)
-                {
-                    using MeshElement = typename ET::Element;
-                    using ShapeFunction2ndOrder = typename ET::ShapeFunction;
-                    _builder[std::type_index(typeid(MeshElement))] =
-                        makeLocalAssemblerBuilder<ShapeFunction2ndOrder>();
-                });
-        }
-    }
-
-    /// Returns data pointer to the newly created local assembler data.
-    ///
-    /// \attention
-    /// The index \c id is not necessarily the mesh item's id. Especially when
-    /// having multiple meshes it will differ from the latter.
-    LADataIntfPtr operator()(std::size_t const id,
-                             MeshLib::Element const& mesh_item,
-                             ConstructorArgs&&... args) const
-    {
-        auto const type_idx = std::type_index(typeid(mesh_item));
-        auto const it = _builder.find(type_idx);
-
-        if (it != _builder.end())
-        {
-            auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
-            return it->second(mesh_item, num_local_dof,
-                              std::forward<ConstructorArgs>(args)...);
-        }
-        OGS_FATAL(
-            "You are trying to build a local assembler for an unknown mesh "
-            "element type ({:s})."
-            " Maybe you have disabled this mesh element type in your build "
-            "configuration, or a mesh element order does not match shape "
-            "function order given in the project file.",
-            type_idx.name());
-    }
-
-private:
-    using LADataBuilder =
-        std::function<LADataIntfPtr(MeshLib::Element const& e,
-                                    std::size_t const local_matrix_size,
-                                    ConstructorArgs&&...)>;
-
-    template <typename ShapeFunction>
-    using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
-        typename ShapeFunction::MeshElement>::IntegrationMethod;
-
-    template <typename ShapeFunction>
-    using LAData =
-        LocalAssemblerData<ShapeFunction, IntegrationMethod<ShapeFunction>,
-                           GlobalDim>;
-
-    /// Generates a function that creates a new LocalAssembler of type
-    /// LAData<ShapeFunction>. Only functions with shape function's dimension
-    /// less or equal to the global dimension are instantiated, e.g.  following
-    /// combinations of shape functions and global dimensions: (Line2, 1),
-    /// (Line2, 2), (Line2, 3), (Hex20, 3) but not (Hex20, 2) or (Hex20, 1).
-    template <typename ShapeFunction>
-    static LADataBuilder makeLocalAssemblerBuilder()
-    {
-        return [](MeshLib::Element const& e,
-                  std::size_t const local_matrix_size,
-                  ConstructorArgs&&... args)
-        {
-            return LADataIntfPtr{new LAData<ShapeFunction>{
-                e, local_matrix_size, std::forward<ConstructorArgs>(args)...}};
-        };
-    }
-
-    /// Mapping of element types to local assembler constructors.
-    std::unordered_map<std::type_index, LADataBuilder> _builder;
-
-    NumLib::LocalToGlobalIndexMap const& _dof_table;
-};
-
-}  // namespace BoundaryConditionAndSourceTerm
-}  // namespace ProcessLib
diff --git a/ProcessLib/HydroMechanics/CreateLocalAssemblers.h b/ProcessLib/HydroMechanics/CreateLocalAssemblers.h
deleted file mode 100644
index ddbb5937fbdc6967100e7d8c915e1914556aae65..0000000000000000000000000000000000000000
--- a/ProcessLib/HydroMechanics/CreateLocalAssemblers.h
+++ /dev/null
@@ -1,84 +0,0 @@
-/**
- * \file
- * \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
-
-#include <vector>
-
-#include "BaseLib/Logging.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "ProcessLib/Utils/LocalDataInitializerHM.h"
-
-namespace ProcessLib
-{
-namespace HydroMechanics
-{
-namespace detail
-{
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    // Shape matrices initializer
-    using LocalDataInitializer =
-        LocalDataInitializerHM<LocalAssemblerInterface,
-                               LocalAssemblerImplementation, GlobalDim,
-                               ExtraCtorArgs...>;
-
-    DBUG("Create local assemblers.");
-    // Populate the vector of local assemblers.
-    local_assemblers.resize(mesh_elements.size());
-
-    LocalDataInitializer initializer(dof_table);
-
-    DBUG("Calling local assembler builder for all mesh elements.");
-    GlobalExecutor::transformDereferenced(
-        initializer, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-
-}  // namespace detail
-
-/*! Creates local assemblers for each element of the given \c mesh.
- *
- * \tparam LocalAssemblerImplementation the individual local assembler type
- * \tparam LocalAssemblerInterface the general local assembler interface
- * \tparam ExtraCtorArgs types of additional constructor arguments.
- *         Those arguments will be passed to the constructor of
- *         \c LocalAssemblerImplementation.
- *
- * The first two template parameters cannot be deduced from the arguments.
- * Therefore they always have to be provided manually.
- */
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    const unsigned /*dimension*/,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    DBUG("Create local assemblers.");
-
-    detail::createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
-        dof_table, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-}  // namespace HydroMechanics
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
index c49de80bcf901875deb788b3002cc530cbcf1334..751c390c02ca3e15e15793f651f2d938e82b085a 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.cpp
@@ -17,8 +17,9 @@
 #include "MeshLib/Elements/Utils.h"
 #include "NumLib/DOF/ComputeSparsityPattern.h"
 #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
-#include "ProcessLib/HydroMechanics/CreateLocalAssemblers.h"
 #include "ProcessLib/Process.h"
+#include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h"
+
 namespace ProcessLib
 {
 namespace HydroMechanics
@@ -170,9 +171,9 @@ void HydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
     MeshLib::Mesh const& mesh,
     unsigned const integration_order)
 {
-    ProcessLib::HydroMechanics::createLocalAssemblers<
-        DisplacementDim, HydroMechanicsLocalAssembler>(
-        mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+    ProcessLib::createLocalAssemblersHM<DisplacementDim,
+                                        HydroMechanicsLocalAssembler>(
+        mesh.getElements(), dof_table, _local_assemblers,
         mesh.isAxiallySymmetric(), integration_order, _process_data);
 
     auto add_secondary_variable = [&](std::string const& name,
diff --git a/ProcessLib/RichardsMechanics/CreateLocalAssemblers.h b/ProcessLib/RichardsMechanics/CreateLocalAssemblers.h
deleted file mode 100644
index a2d563ee488179b6985c3b32a28cabce4b5d0b4d..0000000000000000000000000000000000000000
--- a/ProcessLib/RichardsMechanics/CreateLocalAssemblers.h
+++ /dev/null
@@ -1,84 +0,0 @@
-/**
- * \file
- * \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
-
-#include <vector>
-
-#include "BaseLib/Logging.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "ProcessLib/Utils/LocalDataInitializerHM.h"
-
-namespace ProcessLib
-{
-namespace RichardsMechanics
-{
-namespace detail
-{
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    // Shape matrices initializer
-    using LocalDataInitializer =
-        LocalDataInitializerHM<LocalAssemblerInterface,
-                               LocalAssemblerImplementation, GlobalDim,
-                               ExtraCtorArgs...>;
-
-    DBUG("Create local assemblers.");
-    // Populate the vector of local assemblers.
-    local_assemblers.resize(mesh_elements.size());
-
-    LocalDataInitializer initializer(dof_table);
-
-    DBUG("Calling local assembler builder for all mesh elements.");
-    GlobalExecutor::transformDereferenced(
-        initializer, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-
-}  // namespace detail
-
-/*! Creates local assemblers for each element of the given \c mesh.
- *
- * \tparam LocalAssemblerImplementation the individual local assembler type
- * \tparam LocalAssemblerInterface the general local assembler interface
- * \tparam ExtraCtorArgs types of additional constructor arguments.
- *         Those arguments will be passed to the constructor of
- *         \c LocalAssemblerImplementation.
- *
- * The first two template parameters cannot be deduced from the arguments.
- * Therefore they always have to be provided manually.
- */
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    const unsigned /*dimension*/,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    DBUG("Create local assemblers.");
-
-    detail::createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
-        dof_table, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-}  // namespace RichardsMechanics
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
index cd46be913bfbbea726bb20b34bd1bf7957710494..5fa79b9ad8534a8b2d1adca2591b70b5732d952b 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsProcess.cpp
@@ -15,7 +15,7 @@
 #include "MeshLib/Elements/Utils.h"
 #include "NumLib/DOF/ComputeSparsityPattern.h"
 #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
-#include "ProcessLib/RichardsMechanics/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h"
 #include "RichardsMechanicsFEM.h"
 #include "RichardsMechanicsProcessData.h"
 
@@ -196,9 +196,9 @@ void RichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
 {
     using nlohmann::json;
 
-    ProcessLib::RichardsMechanics::createLocalAssemblers<
-        DisplacementDim, RichardsMechanicsLocalAssembler>(
-        mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+    ProcessLib::createLocalAssemblersHM<DisplacementDim,
+                                        RichardsMechanicsLocalAssembler>(
+        mesh.getElements(), dof_table, _local_assemblers,
         mesh.isAxiallySymmetric(), integration_order, _process_data);
 
     auto add_secondary_variable = [&](std::string const& name,
diff --git a/ProcessLib/SmallDeformation/CreateLocalAssemblers.h b/ProcessLib/SmallDeformation/CreateLocalAssemblers.h
index 42847cd36424906889b951f8de5752b008bc8924..fcf8d3156d0a76a574fdf2dac426c2ece1063338 100644
--- a/ProcessLib/SmallDeformation/CreateLocalAssemblers.h
+++ b/ProcessLib/SmallDeformation/CreateLocalAssemblers.h
@@ -12,8 +12,8 @@
 #include <vector>
 
 #include "BaseLib/Logging.h"
-#include "LocalDataInitializer.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/Utils/LocalAssemblerFactoryForDimGreaterEqualN.h"
 
 namespace ProcessLib
 {
@@ -30,21 +30,19 @@ void createLocalAssemblers(
     std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
     ExtraCtorArgs&&... extra_ctor_args)
 {
-    // Shape matrices initializer
-    using LocalDataInitializer =
-        LocalDataInitializer<LocalAssemblerInterface,
-                             LocalAssemblerImplementation, GlobalDim,
-                             ExtraCtorArgs...>;
+    using LocAsmFactory =
+        ProcessLib::LocalAssemberFactorySD<LocalAssemblerInterface,
+                                           LocalAssemblerImplementation,
+                                           GlobalDim, ExtraCtorArgs...>;
 
     DBUG("Create local assemblers.");
-    // Populate the vector of local assemblers.
-    local_assemblers.resize(mesh_elements.size());
 
-    LocalDataInitializer initializer(dof_table);
+    LocAsmFactory factory(dof_table);
+    local_assemblers.resize(mesh_elements.size());
 
     DBUG("Calling local assembler builder for all mesh elements.");
     GlobalExecutor::transformDereferenced(
-        initializer, mesh_elements, local_assemblers,
+        factory, mesh_elements, local_assemblers,
         std::forward<ExtraCtorArgs>(extra_ctor_args)...);
 }
 
diff --git a/ProcessLib/SmallDeformation/LocalDataInitializer.h b/ProcessLib/SmallDeformation/LocalDataInitializer.h
deleted file mode 100644
index bafa853ccfa4cbf0c4853eba4def23c5f0f4f882..0000000000000000000000000000000000000000
--- a/ProcessLib/SmallDeformation/LocalDataInitializer.h
+++ /dev/null
@@ -1,139 +0,0 @@
-/**
- * \file
- * \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
-
-#include <functional>
-#include <memory>
-#include <type_traits>
-#include <typeindex>
-#include <typeinfo>
-#include <unordered_map>
-
-#include "MeshLib/Elements/Elements.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
-#include "ProcessLib/Utils/EnabledElements.h"
-
-namespace ProcessLib
-{
-/// The LocalDataInitializer is a functor creating a local assembler data with
-/// corresponding to the mesh element type shape functions and calling
-/// initialization of the new local assembler data.
-/// For example for MeshLib::Quad a local assembler data with template argument
-/// NumLib::ShapeQuad4 is created.
-template <typename LocalAssemblerInterface,
-          template <typename, typename, int>
-          class SmallDeformationLocalAssembler,
-          int GlobalDim, typename... ConstructorArgs>
-class LocalDataInitializer final
-{
-    struct IsElementEnabled
-    {
-        template <typename ElementTraits>
-        constexpr bool operator()(ElementTraits*) const
-        {
-            if constexpr (GlobalDim < ElementTraits::ShapeFunction::DIM)
-            {
-                return false;
-            }
-
-            // mechanics processes in OGS are defined in 2D and 3D only
-            return ElementTraits::Element::dimension >= 2;
-        }
-    };
-
-public:
-    using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
-
-    explicit LocalDataInitializer(
-        NumLib::LocalToGlobalIndexMap const& dof_table)
-        : _dof_table(dof_table)
-    {
-        using EnabledElementTraits =
-            decltype(BaseLib::TMP::filter<EnabledElementTraitsLagrange>(
-                std::declval<IsElementEnabled>()));
-
-        BaseLib::TMP::foreach<EnabledElementTraits>(
-            [this]<typename ET>(ET*)
-            {
-                using MeshElement = typename ET::Element;
-                using ShapeFunction = typename ET::ShapeFunction;
-
-                _builder[std::type_index(typeid(MeshElement))] =
-                    makeLocalAssemblerBuilder<ShapeFunction>();
-            });
-    }
-
-    /// Returns data pointer to the newly created local assembler data.
-    ///
-    /// \attention
-    /// The index \c id is not necessarily the mesh item's id. Especially when
-    /// having multiple meshes it will differ from the latter.
-    LADataIntfPtr operator()(std::size_t const id,
-                             MeshLib::Element const& mesh_item,
-                             ConstructorArgs&&... args) const
-    {
-        auto const type_idx = std::type_index(typeid(mesh_item));
-        auto const it = _builder.find(type_idx);
-
-        if (it != _builder.end())
-        {
-            auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
-            return it->second(mesh_item, num_local_dof,
-                              std::forward<ConstructorArgs>(args)...);
-        }
-        OGS_FATAL(
-            "You are trying to build a local assembler for an unknown mesh "
-            "element type ({:s})."
-            " Maybe you have disabled this mesh element type in your build "
-            "configuration, or a mesh element order does not match shape "
-            "function order given in the project file.",
-            type_idx.name());
-    }
-
-private:
-    using LADataBuilder =
-        std::function<LADataIntfPtr(MeshLib::Element const& e,
-                                    std::size_t const local_matrix_size,
-                                    ConstructorArgs&&...)>;
-
-    template <typename ShapeFunction>
-    using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
-        typename ShapeFunction::MeshElement>::IntegrationMethod;
-
-    template <typename ShapeFunction>
-    using LAData = SmallDeformationLocalAssembler<
-        ShapeFunction, IntegrationMethod<ShapeFunction>, GlobalDim>;
-
-    /// Generates a function that creates a new LocalAssembler of type
-    /// LAData<ShapeFunction>. Only functions with shape function's dimension
-    /// less or equal to the global dimension are instantiated, e.g.  following
-    /// combinations of shape functions and global dimensions: (Line2, 1),
-    /// (Line2, 2), (Line2, 3), (Hex20, 3) but not (Hex20, 2) or (Hex20, 1).
-    template <typename ShapeFunction>
-    static LADataBuilder makeLocalAssemblerBuilder()
-    {
-        return [](MeshLib::Element const& e,
-                  std::size_t const local_matrix_size,
-                  ConstructorArgs&&... args)
-        {
-            return LADataIntfPtr{new LAData<ShapeFunction>{
-                e, local_matrix_size, std::forward<ConstructorArgs>(args)...}};
-        };
-    }
-
-    /// Mapping of element types to local assembler constructors.
-    std::unordered_map<std::type_index, LADataBuilder> _builder;
-
-    NumLib::LocalToGlobalIndexMap const& _dof_table;
-};
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/StokesFlow/CreateLocalAssemblers.h b/ProcessLib/StokesFlow/CreateLocalAssemblers.h
deleted file mode 100644
index 887b298963ac335b5ee93d689dc6b775c956cd8c..0000000000000000000000000000000000000000
--- a/ProcessLib/StokesFlow/CreateLocalAssemblers.h
+++ /dev/null
@@ -1,84 +0,0 @@
-/**
- * \file
- * \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
-
-#include <vector>
-
-#include "BaseLib/Logging.h"
-#include "LocalDataInitializer.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-
-namespace ProcessLib
-{
-namespace StokesFlow
-{
-namespace detail
-{
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    // Shape matrices initializer
-    using LocalDataInitializer =
-        LocalDataInitializer<LocalAssemblerInterface,
-                             LocalAssemblerImplementation, GlobalDim,
-                             ExtraCtorArgs...>;
-
-    DBUG("Create local assemblers.");
-    // Populate the vector of local assemblers.
-    local_assemblers.resize(mesh_elements.size());
-
-    LocalDataInitializer initializer(dof_table);
-
-    DBUG("Calling local assembler builder for all mesh elements.");
-    GlobalExecutor::transformDereferenced(
-        initializer, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-
-}  // namespace detail
-
-/*! Creates local assemblers for each element of the given \c mesh.
- *
- * \tparam LocalAssemblerImplementation the individual local assembler type
- * \tparam LocalAssemblerInterface the general local assembler interface
- * \tparam ExtraCtorArgs types of additional constructor arguments.
- *         Those arguments will be passed to the constructor of
- *         \c LocalAssemblerImplementation.
- *
- * The first two template parameters cannot be deduced from the arguments.
- * Therefore they always have to be provided manually.
- */
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    const unsigned /*dimension*/,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    DBUG("Create local assemblers.");
-
-    detail::createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
-        dof_table, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-}  // namespace StokesFlow
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/StokesFlow/LocalDataInitializer.h b/ProcessLib/StokesFlow/LocalDataInitializer.h
deleted file mode 100644
index e317b1698a6053c984934189d24a05493d4388dd..0000000000000000000000000000000000000000
--- a/ProcessLib/StokesFlow/LocalDataInitializer.h
+++ /dev/null
@@ -1,152 +0,0 @@
-/**
- * \file
- * \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
-
-#include <functional>
-#include <memory>
-#include <type_traits>
-#include <typeindex>
-#include <typeinfo>
-#include <unordered_map>
-
-#include "MeshLib/Elements/Elements.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
-#include "ProcessLib/Utils/EnabledElements.h"
-
-namespace ProcessLib::StokesFlow
-{
-/// The LocalDataInitializer is a functor creating a local assembler data with
-/// corresponding to the mesh element type shape functions and calling
-/// initialization of the new local assembler data.
-/// For example for MeshLib::Quad a local assembler data with template argument
-/// NumLib::ShapeQuad4 is created.
-///
-/// \attention This is modified version of the ProcessLib::LocalDataInitializer
-/// class which does not include line or point elements. For the shape functions
-/// of order 2 (used for fluid velocity in StokesFlow) a shape function of order
-/// 1 will be used for the scalar variables.
-template <typename LocalAssemblerInterface,
-          template <typename, typename, typename, int> class LocalAssemblerData,
-          int GlobalDim, typename... ConstructorArgs>
-class LocalDataInitializer final
-{
-    struct IsElementEnabled
-    {
-        template <typename ElementTraits>
-        constexpr bool operator()(ElementTraits*) const
-        {
-            if constexpr (GlobalDim < ElementTraits::ShapeFunction::DIM)
-            {
-                return false;
-            }
-
-            return ElementTraits::Element::dimension >= 2 &&
-                   ElementTraits::ShapeFunction::ORDER == 2;
-        }
-    };
-
-public:
-    using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
-
-    explicit LocalDataInitializer(
-        NumLib::LocalToGlobalIndexMap const& dof_table)
-        : _dof_table(dof_table)
-    {
-        using EnabledElementTraits =
-            decltype(BaseLib::TMP::filter<EnabledElementTraitsLagrange>(
-                std::declval<IsElementEnabled>()));
-
-        BaseLib::TMP::foreach<EnabledElementTraits>(
-            [this]<typename ET>(ET*)
-            {
-                using MeshElement = typename ET::Element;
-                using ShapeFunction = typename ET::ShapeFunction;
-                using LowerOrderShapeFunction =
-                    typename ET::LowerOrderShapeFunction;
-
-                _builder[std::type_index(typeid(MeshElement))] =
-                    makeLocalAssemblerBuilder<ShapeFunction,
-                                              LowerOrderShapeFunction>();
-            });
-    }
-
-    /// Returns data pointer to the newly created local assembler data.
-    ///
-    /// \attention
-    /// The index \c id is not necessarily the mesh item's id. Especially when
-    /// having multiple meshes it will differ from the latter.
-    LADataIntfPtr operator()(std::size_t const id,
-                             MeshLib::Element const& mesh_item,
-                             ConstructorArgs&&... args) const
-    {
-        auto const type_idx = std::type_index(typeid(mesh_item));
-        auto const it = _builder.find(type_idx);
-
-        if (it != _builder.end())
-        {
-            auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
-            return it->second(mesh_item, num_local_dof,
-                              std::forward<ConstructorArgs>(args)...);
-        }
-        OGS_FATAL(
-            "You are trying to build a local assembler for an unknown mesh "
-            "element type ({:s})."
-            " Maybe you have disabled this mesh element type in your build "
-            "configuration, or a mesh element order does not match shape "
-            "function order given in the project file.",
-            type_idx.name());
-    }
-
-private:
-    using LADataBuilder =
-        std::function<LADataIntfPtr(MeshLib::Element const& e,
-                                    std::size_t const local_matrix_size,
-                                    ConstructorArgs&&...)>;
-
-    template <typename ShapeFunction>
-    using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
-        typename ShapeFunction::MeshElement>::IntegrationMethod;
-
-    template <typename ShapeFunctionDisplacement,
-              typename ShapeFunctionPressure>
-    using LAData =
-        LocalAssemblerData<ShapeFunctionDisplacement, ShapeFunctionPressure,
-                           IntegrationMethod<ShapeFunctionDisplacement>,
-                           GlobalDim>;
-
-    /// Generates a function that creates a new LocalAssembler of type
-    /// LAData<ShapeFunctionDisplacement, ShapeFunctionPressure>. Only functions
-    /// with shape function's dimension less or equal to the global dimension
-    /// are instantiated, e.g.  following combinations of shape functions and
-    /// global dimensions: (Line2, 1), (Line2, 2), (Line2, 3), (Hex20, 3) but
-    /// not (Hex20, 2) or (Hex20, 1).
-    template <typename ShapeFunction, typename LowerOrderShapeFunction>
-    static LADataBuilder makeLocalAssemblerBuilder()
-    {
-        return [](MeshLib::Element const& e,
-                  std::size_t const local_matrix_size,
-                  ConstructorArgs&&... args)
-        {
-            return LADataIntfPtr{
-                new LAData<ShapeFunction, LowerOrderShapeFunction>{
-                    e, local_matrix_size,
-                    std::forward<ConstructorArgs>(args)...}};
-        };
-    }
-
-    /// Mapping of element types to local assembler constructors.
-    std::unordered_map<std::type_index, LADataBuilder> _builder;
-
-    NumLib::LocalToGlobalIndexMap const& _dof_table;
-};
-
-}  // namespace ProcessLib::StokesFlow
diff --git a/ProcessLib/StokesFlow/StokesFlowProcess.cpp b/ProcessLib/StokesFlow/StokesFlowProcess.cpp
index 49b65f9474d714d79d61041900280618de8e0eac..21060f27bb5b24c114334315bcb957b2f8986a64 100644
--- a/ProcessLib/StokesFlow/StokesFlowProcess.cpp
+++ b/ProcessLib/StokesFlow/StokesFlowProcess.cpp
@@ -12,8 +12,8 @@
 
 #include <cassert>
 
-#include "CreateLocalAssemblers.h"
 #include "MeshLib/Elements/Utils.h"
+#include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h"
 
 namespace ProcessLib
 {
@@ -88,12 +88,8 @@ void StokesFlowProcess<GlobalDim>::initializeConcreteProcess(
     MeshLib::Mesh const& mesh,
     unsigned const integration_order)
 {
-    // Todo: may move LocalDataInitializer.h and CreateLocalAssemblers.h which
-    // are identical to those in such processes as HydroMechanics,
-    // RichardsMechanics, and etc, into a common place.
-    ProcessLib::StokesFlow::createLocalAssemblers<GlobalDim,
-                                                  LocalAssemblerData>(
-        mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+    ProcessLib::createLocalAssemblersStokes<GlobalDim, LocalAssemblerData>(
+        mesh.getElements(), dof_table, _local_assemblers,
         mesh.isAxiallySymmetric(), integration_order, _process_data);
 
     _process_data.pressure_interpolated =
diff --git a/ProcessLib/TH2M/CreateLocalAssemblers.h b/ProcessLib/TH2M/CreateLocalAssemblers.h
deleted file mode 100644
index 508d863d2dc358a6f88d4f61f7bbb5f9a3d8559e..0000000000000000000000000000000000000000
--- a/ProcessLib/TH2M/CreateLocalAssemblers.h
+++ /dev/null
@@ -1,84 +0,0 @@
-/**
- * \file
- * \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
-
-#include <vector>
-
-#include "BaseLib/Logging.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "ProcessLib/Utils/LocalDataInitializerHM.h"
-
-namespace ProcessLib
-{
-namespace TH2M
-{
-namespace detail
-{
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    // Shape matrices initializer
-    using LocalDataInitializer =
-        LocalDataInitializerHM<LocalAssemblerInterface,
-                               LocalAssemblerImplementation, GlobalDim,
-                               ExtraCtorArgs...>;
-
-    DBUG("Create local assemblers.");
-    // Populate the vector of local assemblers.
-    local_assemblers.resize(mesh_elements.size());
-
-    LocalDataInitializer initializer(dof_table);
-
-    DBUG("Calling local assembler builder for all mesh elements.");
-    GlobalExecutor::transformDereferenced(
-        initializer, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-
-}  // namespace detail
-
-/*! Creates local assemblers for each element of the given \c mesh.
- *
- * \tparam LocalAssemblerImplementation the individual local assembler type
- * \tparam LocalAssemblerInterface the general local assembler interface
- * \tparam ExtraCtorArgs types of additional constructor arguments.
- *         Those arguments will be passed to the constructor of
- *         \c LocalAssemblerImplementation.
- *
- * The first two template parameters cannot be deduced from the arguments.
- * Therefore they always have to be provided manually.
- */
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    const unsigned /*dimension*/,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    DBUG("Create local assemblers.");
-
-    detail::createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
-        dof_table, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-}  // namespace TH2M
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/TH2M/TH2MProcess.cpp b/ProcessLib/TH2M/TH2MProcess.cpp
index c174567ce32cd6c640e9220dfb2bf366358aea49..8d74851ae4e6f984666782320366d49ccd17dcc0 100644
--- a/ProcessLib/TH2M/TH2MProcess.cpp
+++ b/ProcessLib/TH2M/TH2MProcess.cpp
@@ -15,7 +15,7 @@
 #include "MeshLib/Elements/Utils.h"
 #include "NumLib/DOF/ComputeSparsityPattern.h"
 #include "ProcessLib/Process.h"
-#include "ProcessLib/TH2M/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h"
 #include "TH2MFEM.h"
 #include "TH2MProcessData.h"
 
@@ -158,9 +158,8 @@ void TH2MProcess<DisplacementDim>::initializeConcreteProcess(
     MeshLib::Mesh const& mesh,
     unsigned const integration_order)
 {
-    ProcessLib::TH2M::createLocalAssemblers<DisplacementDim,
-                                            TH2MLocalAssembler>(
-        mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+    ProcessLib::createLocalAssemblersHM<DisplacementDim, TH2MLocalAssembler>(
+        mesh.getElements(), dof_table, _local_assemblers,
         mesh.isAxiallySymmetric(), integration_order, _process_data);
 
     auto add_secondary_variable = [&](std::string const& name,
diff --git a/ProcessLib/ThermoHydroMechanics/CreateLocalAssemblers.h b/ProcessLib/ThermoHydroMechanics/CreateLocalAssemblers.h
deleted file mode 100644
index d512d6de635908378e4b3b26743485020ce45c19..0000000000000000000000000000000000000000
--- a/ProcessLib/ThermoHydroMechanics/CreateLocalAssemblers.h
+++ /dev/null
@@ -1,84 +0,0 @@
-/**
- * \file
- * \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
-
-#include <vector>
-
-#include "BaseLib/Logging.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "ProcessLib/Utils/LocalDataInitializerHM.h"
-
-namespace ProcessLib
-{
-namespace ThermoHydroMechanics
-{
-namespace detail
-{
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    // Shape matrices initializer
-    using LocalDataInitializer =
-        LocalDataInitializerHM<LocalAssemblerInterface,
-                               LocalAssemblerImplementation, GlobalDim,
-                               ExtraCtorArgs...>;
-
-    DBUG("Create local assemblers.");
-    // Populate the vector of local assemblers.
-    local_assemblers.resize(mesh_elements.size());
-
-    LocalDataInitializer initializer(dof_table);
-
-    DBUG("Calling local assembler builder for all mesh elements.");
-    GlobalExecutor::transformDereferenced(
-        initializer, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-
-}  // namespace detail
-
-/*! Creates local assemblers for each element of the given \c mesh.
- *
- * \tparam LocalAssemblerImplementation the individual local assembler type
- * \tparam LocalAssemblerInterface the general local assembler interface
- * \tparam ExtraCtorArgs types of additional constructor arguments.
- *         Those arguments will be passed to the constructor of
- *         \c LocalAssemblerImplementation.
- *
- * The first two template parameters cannot be deduced from the arguments.
- * Therefore they always have to be provided manually.
- */
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    const unsigned /*dimension*/,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    DBUG("Create local assemblers.");
-
-    detail::createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
-        dof_table, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-}  // namespace ThermoHydroMechanics
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
index cc7a4869250c3ea7cdc389d59084901f9e2a724e..d1fd8f1fded488be7eb507e3788ef7d1d8bddd68 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsProcess.cpp
@@ -15,7 +15,7 @@
 #include "MeshLib/Elements/Utils.h"
 #include "NumLib/DOF/ComputeSparsityPattern.h"
 #include "ProcessLib/Process.h"
-#include "ProcessLib/ThermoHydroMechanics/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h"
 #include "ThermoHydroMechanicsFEM.h"
 #include "ThermoHydroMechanicsProcessData.h"
 
@@ -176,9 +176,9 @@ void ThermoHydroMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
     MeshLib::Mesh const& mesh,
     unsigned const integration_order)
 {
-    ProcessLib::ThermoHydroMechanics::createLocalAssemblers<
-        DisplacementDim, ThermoHydroMechanicsLocalAssembler>(
-        mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+    ProcessLib::createLocalAssemblersHM<DisplacementDim,
+                                        ThermoHydroMechanicsLocalAssembler>(
+        mesh.getElements(), dof_table, _local_assemblers,
         mesh.isAxiallySymmetric(), integration_order, _process_data);
 
     _secondary_variables.addSecondaryVariable(
diff --git a/ProcessLib/ThermoRichardsMechanics/CreateLocalAssemblers.h b/ProcessLib/ThermoRichardsMechanics/CreateLocalAssemblers.h
deleted file mode 100644
index ab62ac3bd7f503155c65b3bb1d4190691b81be58..0000000000000000000000000000000000000000
--- a/ProcessLib/ThermoRichardsMechanics/CreateLocalAssemblers.h
+++ /dev/null
@@ -1,80 +0,0 @@
-/**
- * \file
- * \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
-
-#include <vector>
-
-#include "BaseLib/Logging.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "ProcessLib/Utils/LocalDataInitializerHM.h"
-
-namespace ProcessLib::ThermoRichardsMechanics
-{
-namespace detail
-{
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    // Shape matrices initializer
-    using LocalDataInitializer =
-        LocalDataInitializerHM<LocalAssemblerInterface,
-                               LocalAssemblerImplementation, GlobalDim,
-                               ExtraCtorArgs...>;
-
-    DBUG("Create local assemblers.");
-    // Populate the vector of local assemblers.
-    local_assemblers.resize(mesh_elements.size());
-
-    LocalDataInitializer initializer(dof_table);
-
-    DBUG("Calling local assembler builder for all mesh elements.");
-    GlobalExecutor::transformDereferenced(
-        initializer, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-
-}  // namespace detail
-
-/*! Creates local assemblers for each element of the given \c mesh.
- *
- * \tparam LocalAssemblerImplementation the individual local assembler type
- * \tparam LocalAssemblerInterface the general local assembler interface
- * \tparam ExtraCtorArgs types of additional constructor arguments.
- *         Those arguments will be passed to the constructor of
- *         \c LocalAssemblerImplementation.
- *
- * The first two template parameters cannot be deduced from the arguments.
- * Therefore they always have to be provided manually.
- */
-template <int GlobalDim,
-          template <typename, typename, typename, int>
-          class LocalAssemblerImplementation,
-          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
-void createLocalAssemblers(
-    const unsigned /*dimension*/,
-    std::vector<MeshLib::Element*> const& mesh_elements,
-    NumLib::LocalToGlobalIndexMap const& dof_table,
-    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
-    ExtraCtorArgs&&... extra_ctor_args)
-{
-    DBUG("Create local assemblers.");
-
-    detail::createLocalAssemblers<GlobalDim, LocalAssemblerImplementation>(
-        dof_table, mesh_elements, local_assemblers,
-        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
-}
-}  // namespace ProcessLib::ThermoRichardsMechanics
diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
index cedca58f259b504a7b8479832b81c472e3f51255..1bc652a2c0298702745e37c5d97499df9ef099d5 100644
--- a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
+++ b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsProcess.cpp
@@ -13,11 +13,11 @@
 #include <cassert>
 
 #include "BaseLib/Error.h"
-#include "CreateLocalAssemblers.h"
 #include "MeshLib/Elements/Utils.h"
 #include "NumLib/DOF/ComputeSparsityPattern.h"
 #include "ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h"
 #include "ProcessLib/Process.h"
+#include "ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h"
 #include "ThermoRichardsMechanicsFEM.h"
 
 namespace ProcessLib
@@ -154,9 +154,9 @@ void ThermoRichardsMechanicsProcess<DisplacementDim>::initializeConcreteProcess(
     MeshLib::Mesh const& mesh,
     unsigned const integration_order)
 {
-    createLocalAssemblers<DisplacementDim,
-                          ThermoRichardsMechanicsLocalAssembler>(
-        mesh.getDimension(), mesh.getElements(), dof_table, local_assemblers_,
+    createLocalAssemblersHM<DisplacementDim,
+                            ThermoRichardsMechanicsLocalAssembler>(
+        mesh.getElements(), dof_table, local_assemblers_,
         mesh.isAxiallySymmetric(), integration_order, process_data_);
 
     auto add_secondary_variable = [&](std::string const& name,
diff --git a/ProcessLib/Utils/CreateLocalAssemblers.h b/ProcessLib/Utils/CreateLocalAssemblers.h
index 87288a33994966177754b167b8026f497167f506..37f964ac58454fcdc1e1829d3fb73c989dd08aae 100644
--- a/ProcessLib/Utils/CreateLocalAssemblers.h
+++ b/ProcessLib/Utils/CreateLocalAssemblers.h
@@ -12,7 +12,7 @@
 #include <vector>
 
 #include "BaseLib/Logging.h"
-#include "LocalDataInitializer.h"
+#include "LocalAssemblerFactoryForDimGreaterEqualN.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
 
 namespace ProcessLib
@@ -29,21 +29,19 @@ void createLocalAssemblers(
     static_assert(
         GlobalDim == 1 || GlobalDim == 2 || GlobalDim == 3,
         "Meshes with dimension greater than three are not supported.");
-    // Shape matrices initializer
-    using LocalDataInitializer =
-        LocalDataInitializer<LocalAssemblerInterface,
-                             LocalAssemblerImplementation, GlobalDim,
-                             ExtraCtorArgs...>;
+
+    using LocAsmFactory = LocalAssemberFactory<LocalAssemblerInterface,
+                                               LocalAssemblerImplementation,
+                                               GlobalDim, ExtraCtorArgs...>;
 
     DBUG("Create local assemblers.");
-    // Populate the vector of local assemblers.
-    local_assemblers.resize(mesh_elements.size());
 
-    LocalDataInitializer initializer(dof_table);
+    LocAsmFactory factory(dof_table);
+    local_assemblers.resize(mesh_elements.size());
 
     DBUG("Calling local assembler builder for all mesh elements.");
     GlobalExecutor::transformDereferenced(
-        initializer, mesh_elements, local_assemblers,
+        factory, mesh_elements, local_assemblers,
         std::forward<ExtraCtorArgs>(extra_ctor_args)...);
 }
 
diff --git a/ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h b/ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h
new file mode 100644
index 0000000000000000000000000000000000000000..99641df7ede28a74b166a9cacfc441701028cba9
--- /dev/null
+++ b/ProcessLib/Utils/CreateLocalAssemblersTaylorHood.h
@@ -0,0 +1,104 @@
+/**
+ * \file
+ * \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
+
+#include <vector>
+
+#include "BaseLib/Logging.h"
+#include "LocalAssemblerFactoryTaylorHoodForOrderGreaterEqualN.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+namespace ProcessLib
+{
+namespace detail
+{
+/*! Creates local assemblers for each element of the given \c mesh.
+ *
+ * \tparam LocalAssemblerFactory the factory that will instantiate the local
+ *         assemblers
+ * \tparam LocalAssemblerImplementation the individual local assembler type
+ * \tparam LocalAssemblerInterface the general local assembler interface
+ * \tparam ExtraCtorArgs types of additional constructor arguments.
+ *         Those arguments will be passed to the constructor of
+ *         \c LocalAssemblerImplementation.
+ */
+template <template <typename /*LocalAssemblerInterface*/,
+                    template <typename /* shp fct */,
+                              typename /* lower order shp fct */,
+                              typename /* int meth */, int /* global dim */>
+                    class /*LocalAssemblerImplementation*/,
+                    int /* global dim */, typename... /*ConstructorArgs*/>
+          class LocalAssemblerFactory,
+          int GlobalDim,
+          template <typename /* shp fct */, typename /* lower order shp fct */,
+                    typename /* int meth */, int /* global dim */>
+          class LocalAssemblerImplementation,
+          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblersTaylorHood(
+    std::vector<MeshLib::Element*> const& mesh_elements,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+    ExtraCtorArgs&&... extra_ctor_args)
+{
+    using LocAsmFac = LocalAssemblerFactory<LocalAssemblerInterface,
+                                            LocalAssemblerImplementation,
+                                            GlobalDim, ExtraCtorArgs...>;
+
+    DBUG("Create local assemblers.");
+
+    LocAsmFac factory(dof_table);
+    local_assemblers.resize(mesh_elements.size());
+
+    DBUG("Calling local assembler builder for all mesh elements.");
+    GlobalExecutor::transformDereferenced(
+        factory, mesh_elements, local_assemblers,
+        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+
+}  // namespace detail
+
+template <int GlobalDim,
+          template <typename /* shp fct */, typename /* lower order shp fct */,
+                    typename /* int meth */, int /* global dim */>
+          class LocalAssemblerImplementation,
+          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblersHM(
+    std::vector<MeshLib::Element*> const& mesh_elements,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+    ExtraCtorArgs&&... extra_ctor_args)
+{
+    detail::createLocalAssemblersTaylorHood<LocalAssemblerFactoryHM, GlobalDim,
+                                            LocalAssemblerImplementation,
+                                            LocalAssemblerInterface>(
+        mesh_elements, dof_table, local_assemblers,
+        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+
+template <int GlobalDim,
+          template <typename /* shp fct */, typename /* lower order shp fct */,
+                    typename /* int meth */, int /* global dim */>
+          class LocalAssemblerImplementation,
+          typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+void createLocalAssemblersStokes(
+    std::vector<MeshLib::Element*> const& mesh_elements,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
+    ExtraCtorArgs&&... extra_ctor_args)
+{
+    detail::createLocalAssemblersTaylorHood<
+        LocalAssemblerFactoryStokes, GlobalDim, LocalAssemblerImplementation,
+        LocalAssemblerInterface>(
+        mesh_elements, dof_table, local_assemblers,
+        std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+}
+
+}  // namespace ProcessLib
diff --git a/ProcessLib/Utils/GenericLocalAssemblerFactory.h b/ProcessLib/Utils/GenericLocalAssemblerFactory.h
new file mode 100644
index 0000000000000000000000000000000000000000..30ddf38414f8bc1ce243565e4ccf10a0f5c4cef5
--- /dev/null
+++ b/ProcessLib/Utils/GenericLocalAssemblerFactory.h
@@ -0,0 +1,113 @@
+/**
+ * \file
+ * \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
+
+#include <functional>
+#include <memory>
+#include <typeindex>
+#include <unordered_map>
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
+
+namespace ProcessLib
+{
+/// A functor creating a local assembler with shape functions corresponding to
+/// the mesh element type.
+template <typename LocalAssemblerInterface, typename... ConstructorArgs>
+struct GenericLocalAssemblerFactory
+{
+    using LocAsmIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
+    using LocAsmBuilder =
+        std::function<LocAsmIntfPtr(MeshLib::Element const& e,
+                                    std::size_t const local_matrix_size,
+                                    ConstructorArgs&&...)>;
+
+protected:  // only allow instances of subclasses
+    explicit GenericLocalAssemblerFactory(
+        NumLib::LocalToGlobalIndexMap const& dof_table)
+        : _dof_table(dof_table)
+    {
+    }
+
+public:
+    /// Returns the newly created local assembler.
+    ///
+    /// \attention
+    /// The index \c id is not necessarily the mesh item's id. Especially when
+    /// having multiple meshes it will differ from the latter.
+    LocAsmIntfPtr operator()(std::size_t const id,
+                             MeshLib::Element const& mesh_item,
+                             ConstructorArgs&&... args) const
+    {
+        auto const type_idx = std::type_index(typeid(mesh_item));
+        auto const it = _builders.find(type_idx);
+
+        if (it != _builders.end())
+        {
+            auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
+            return it->second(mesh_item, num_local_dof,
+                              std::forward<ConstructorArgs>(args)...);
+        }
+        OGS_FATAL(
+            "You are trying to build a local assembler for an unknown mesh "
+            "element type ({:s})."
+            " Maybe you have disabled this mesh element type in your build "
+            "configuration, or a mesh element order does not match shape "
+            "function order given in the project file.",
+            type_idx.name());
+    }
+
+private:
+    NumLib::LocalToGlobalIndexMap const& _dof_table;
+
+protected:
+    /// Mapping of element types to local assembler builders.
+    std::unordered_map<std::type_index, LocAsmBuilder> _builders;
+};
+
+template <typename ShapeFunction, typename LocalAssemblerInterface,
+          template <typename /* shp fct */, typename /* int meth */,
+                    int /* global dim */>
+          class LocalAssemblerImplementation,
+          int GlobalDim, typename... ConstructorArgs>
+class LocalAssemblerBuilderFactory
+{
+    using GLAF = GenericLocalAssemblerFactory<LocalAssemblerInterface,
+                                              ConstructorArgs...>;
+    using LocAsmIntfPtr = typename GLAF::LocAsmIntfPtr;
+    using LocAsmBuilder = typename GLAF::LocAsmBuilder;
+
+    using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
+        typename ShapeFunction::MeshElement>::IntegrationMethod;
+
+    using LocAsmImpl =
+        LocalAssemblerImplementation<ShapeFunction, IntegrationMethod,
+                                     GlobalDim>;
+
+    LocalAssemblerBuilderFactory() = delete;
+
+public:
+    /// Generates a function that creates a new local assembler of type
+    /// \c LocAsmImpl.
+    static LocAsmBuilder create()
+    {
+        return [](MeshLib::Element const& e,
+                  std::size_t const local_matrix_size,
+                  ConstructorArgs&&... args)
+        {
+            return std::make_unique<LocAsmImpl>(
+                e, local_matrix_size, std::forward<ConstructorArgs>(args)...);
+        };
+    }
+};
+
+}  // namespace ProcessLib
diff --git a/ProcessLib/Utils/LocalAssemblerFactoryForDimGreaterEqualN.h b/ProcessLib/Utils/LocalAssemblerFactoryForDimGreaterEqualN.h
new file mode 100644
index 0000000000000000000000000000000000000000..3778e37b51726bedc33fd52e9deeda34446082ad
--- /dev/null
+++ b/ProcessLib/Utils/LocalAssemblerFactoryForDimGreaterEqualN.h
@@ -0,0 +1,101 @@
+/**
+ * \file
+ * \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
+
+#include "EnabledElements.h"
+#include "GenericLocalAssemblerFactory.h"
+
+namespace ProcessLib
+{
+/**
+ * Creates local assemblers for mesh elements with dimension greater than or
+ * equal to \c MinElementDim.
+ *
+ * This local assembler factory supports a single type of shape functions, i.e.,
+ * all primary variables are discretized with the same shape function.
+ */
+template <int MinElementDim,
+          typename LocalAssemblerInterface,
+          template <typename, typename, int>
+          class LocalAssemblerImplementation,
+          int GlobalDim,
+          typename... ConstructorArgs>
+class LocalAssemblerFactoryForDimGreaterEqualN final
+    : public GenericLocalAssemblerFactory<LocalAssemblerInterface,
+                                          ConstructorArgs...>
+{
+    using Base = GenericLocalAssemblerFactory<LocalAssemblerInterface,
+                                              ConstructorArgs...>;
+
+    struct IsElementEnabled
+    {
+        template <typename ElementTraits>
+        constexpr bool operator()(ElementTraits*) const
+        {
+            if constexpr (GlobalDim < ElementTraits::ShapeFunction::DIM)
+            {
+                return false;
+            }
+
+            return ElementTraits::Element::dimension >= MinElementDim;
+        }
+    };
+
+public:
+    explicit LocalAssemblerFactoryForDimGreaterEqualN(
+        NumLib::LocalToGlobalIndexMap const& dof_table)
+        : Base{dof_table}
+    {
+        using EnabledElementTraits =
+            decltype(BaseLib::TMP::filter<EnabledElementTraitsLagrange>(
+                std::declval<IsElementEnabled>()));
+
+        BaseLib::TMP::foreach<EnabledElementTraits>(
+            [this]<typename ET>(ET*)
+            {
+                using MeshElement = typename ET::Element;
+                using ShapeFunction = typename ET::ShapeFunction;
+                Base::_builders[std::type_index(typeid(MeshElement))] =
+                    LocalAssemblerBuilderFactory<ShapeFunction,
+                                                 LocalAssemblerInterface,
+                                                 LocalAssemblerImplementation,
+                                                 GlobalDim,
+                                                 ConstructorArgs...>::create();
+            });
+    }
+};
+
+/// By default processes in OGS are defined in 1D, 2D and 3D.
+template <typename LocalAssemblerInterface,
+          template <typename, typename, int>
+          class LocalAssemblerImplementation,
+          int GlobalDim,
+          typename... ConstructorArgs>
+using LocalAssemberFactory =
+    LocalAssemblerFactoryForDimGreaterEqualN<1,
+                                             LocalAssemblerInterface,
+                                             LocalAssemblerImplementation,
+                                             GlobalDim,
+                                             ConstructorArgs...>;
+
+/// Mechanics processes in OGS are defined in 2D and 3D only.
+template <typename LocalAssemblerInterface,
+          template <typename, typename, int>
+          class LocalAssemblerImplementation,
+          int GlobalDim,
+          typename... ConstructorArgs>
+using LocalAssemberFactorySD =
+    LocalAssemblerFactoryForDimGreaterEqualN<2,
+                                             LocalAssemblerInterface,
+                                             LocalAssemblerImplementation,
+                                             GlobalDim,
+                                             ConstructorArgs...>;
+}  // namespace ProcessLib
diff --git a/ProcessLib/Utils/LocalAssemblerFactoryTaylorHoodForOrderGreaterEqualN.h b/ProcessLib/Utils/LocalAssemblerFactoryTaylorHoodForOrderGreaterEqualN.h
new file mode 100644
index 0000000000000000000000000000000000000000..d3ca34e531b88c67307a72128dccb702b737316b
--- /dev/null
+++ b/ProcessLib/Utils/LocalAssemblerFactoryTaylorHoodForOrderGreaterEqualN.h
@@ -0,0 +1,174 @@
+/**
+ * \file
+ * \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
+
+#include "EnabledElements.h"
+#include "GenericLocalAssemblerFactory.h"
+
+namespace ProcessLib
+{
+template <typename ShapeFunction,
+          typename LowerOrderShapeFunction,
+          typename LocalAssemblerInterface,
+          template <typename /* shp fct */,
+                    typename /* lower order shp fct */,
+                    typename /* int meth */,
+                    int /* global dim */>
+          class LocalAssemblerImplementation,
+          int GlobalDim,
+          typename... ConstructorArgs>
+class LocalAssemblerBuilderFactoryTaylorHood
+{
+    using GLAF = GenericLocalAssemblerFactory<LocalAssemblerInterface,
+                                              ConstructorArgs...>;
+    using LocAsmIntfPtr = typename GLAF::LocAsmIntfPtr;
+    using LocAsmBuilder = typename GLAF::LocAsmBuilder;
+
+    using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
+        typename ShapeFunction::MeshElement>::IntegrationMethod;
+
+    using LocAsmImpl = LocalAssemblerImplementation<ShapeFunction,
+                                                    LowerOrderShapeFunction,
+                                                    IntegrationMethod,
+                                                    GlobalDim>;
+
+    LocalAssemblerBuilderFactoryTaylorHood() = delete;
+
+public:
+    /// Generates a function that creates a new local assembler of type
+    /// \c LocAsmImpl.
+    static LocAsmBuilder create()
+    {
+        return [](MeshLib::Element const& e,
+                  std::size_t const local_matrix_size,
+                  ConstructorArgs&&... args)
+        {
+            return std::make_unique<LocAsmImpl>(
+                e, local_matrix_size, std::forward<ConstructorArgs>(args)...);
+        };
+    }
+};
+
+/**
+ * Local assembler factory for Taylor-Hood elements.
+ *
+ * Elements/shape functions must be of order greater than or equal to \c
+ * MinShapeFctOrder.
+ *
+ * If \c MinShapeFctOrder is 1, local assemblers are instantiated also for
+ * linear mesh elements. In this case we don't have Taylor-Hood elements for
+ * linear mesh elements. Instead, on linear mesh elements all shape functions
+ * will have the same order (namely 1).
+ */
+template <int MinShapeFctOrder,
+          typename LocalAssemblerInterface,
+          template <typename /* shp fct */,
+                    typename /* lower order shp fct */,
+                    typename /* int meth */,
+                    int /* global dim */>
+          class LocalAssemblerImplementation,
+          int GlobalDim,
+          typename... ConstructorArgs>
+class LocalAssemblerFactoryTaylorHoodForOrderGreaterEqualN final
+    : public ProcessLib::GenericLocalAssemblerFactory<LocalAssemblerInterface,
+                                                      ConstructorArgs...>
+{
+    using Base =
+        ProcessLib::GenericLocalAssemblerFactory<LocalAssemblerInterface,
+                                                 ConstructorArgs...>;
+
+    template <typename ShapeFunction, typename LowerOrderShapeFunction>
+    using LocAsmBuilderFactory =
+        LocalAssemblerBuilderFactoryTaylorHood<ShapeFunction,
+                                               LowerOrderShapeFunction,
+                                               LocalAssemblerInterface,
+                                               LocalAssemblerImplementation,
+                                               GlobalDim,
+                                               ConstructorArgs...>;
+
+    struct IsElementEnabled
+    {
+        template <typename ElementTraits>
+        constexpr bool operator()(ElementTraits*) const
+        {
+            if constexpr (GlobalDim < ElementTraits::ShapeFunction::DIM)
+            {
+                return false;
+            }
+
+            if constexpr (ElementTraits::Element::dimension < 2)
+            {
+                // in OGS Taylor-Hood elements are available in 2D and 3D only
+                return false;
+            }
+
+            return ElementTraits::ShapeFunction::ORDER >= MinShapeFctOrder;
+        }
+    };
+
+public:
+    explicit LocalAssemblerFactoryTaylorHoodForOrderGreaterEqualN(
+        NumLib::LocalToGlobalIndexMap const& dof_table)
+        : Base(dof_table)
+    {
+        using EnabledElementTraits =
+            decltype(BaseLib::TMP::filter<EnabledElementTraitsLagrange>(
+                std::declval<IsElementEnabled>()));
+
+        BaseLib::TMP::foreach<EnabledElementTraits>(
+            [this]<typename ET>(ET*)
+            {
+                using MeshElement = typename ET::Element;
+                using ShapeFunction = typename ET::ShapeFunction;
+                using LowerOrderShapeFunction =
+                    typename ET::LowerOrderShapeFunction;
+
+                Base::_builders[std::type_index(typeid(MeshElement))] =
+                    LocalAssemblerBuilderFactoryTaylorHood<
+                        ShapeFunction,
+                        LowerOrderShapeFunction,
+                        LocalAssemblerInterface,
+                        LocalAssemblerImplementation,
+                        GlobalDim,
+                        ConstructorArgs...>::create();
+            });
+    }
+};
+
+/// HM processes in OGS are defined for linear and higher order elements.
+template <typename LocalAssemblerInterface,
+          template <typename, typename, typename, int>
+          class LocalAssemblerImplementation,
+          int GlobalDim,
+          typename... ConstructorArgs>
+using LocalAssemblerFactoryHM =
+    LocalAssemblerFactoryTaylorHoodForOrderGreaterEqualN<
+        1,
+        LocalAssemblerInterface,
+        LocalAssemblerImplementation,
+        GlobalDim,
+        ConstructorArgs...>;
+
+/// Stokes flow in OGS is defined for higher order elements only.
+template <typename LocalAssemblerInterface,
+          template <typename, typename, typename, int>
+          class LocalAssemblerImplementation,
+          int GlobalDim,
+          typename... ConstructorArgs>
+using LocalAssemblerFactoryStokes =
+    LocalAssemblerFactoryTaylorHoodForOrderGreaterEqualN<
+        2,
+        LocalAssemblerInterface,
+        LocalAssemblerImplementation,
+        GlobalDim,
+        ConstructorArgs...>;
+
+}  // namespace ProcessLib
diff --git a/ProcessLib/Utils/LocalDataInitializer.h b/ProcessLib/Utils/LocalDataInitializer.h
deleted file mode 100644
index 53ae6944e8a91468887976aa997dd97367d88b9a..0000000000000000000000000000000000000000
--- a/ProcessLib/Utils/LocalDataInitializer.h
+++ /dev/null
@@ -1,137 +0,0 @@
-/**
- * \file
- * \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
-
-#include <functional>
-#include <memory>
-#include <type_traits>
-#include <typeindex>
-#include <typeinfo>
-#include <unordered_map>
-
-#include "EnabledElements.h"
-#include "MeshLib/Elements/Elements.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
-
-namespace ProcessLib
-{
-/// The LocalDataInitializer is a functor creating a local assembler data with
-/// corresponding to the mesh element type shape functions and calling
-/// initialization of the new local assembler data.
-/// For example for MeshLib::Quad a local assembler data with template argument
-/// NumLib::ShapeQuad4 is created.
-template <typename LocalAssemblerInterface,
-          template <typename, typename, int> class LocalAssemblerData,
-          int GlobalDim, typename... ConstructorArgs>
-class LocalDataInitializer final
-{
-    struct IsElementEnabled
-    {
-        template <typename ElementTraits>
-        constexpr bool operator()(ElementTraits*) const
-        {
-            if constexpr (GlobalDim < ElementTraits::ShapeFunction::DIM)
-            {
-                return false;
-            }
-
-            // exclude 0D elements
-            return ElementTraits::Element::dimension >= 1;
-        }
-    };
-
-public:
-    using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
-
-    explicit LocalDataInitializer(NumLib::LocalToGlobalIndexMap const& dof_table)
-        : _dof_table(dof_table)
-    {
-        using EnabledElementTraits =
-            decltype(BaseLib::TMP::filter<EnabledElementTraitsLagrange>(
-                std::declval<IsElementEnabled>()));
-
-        BaseLib::TMP::foreach<EnabledElementTraits>(
-            [this]<typename ET>(ET*)
-            {
-                using MeshElement = typename ET::Element;
-                using ShapeFunction = typename ET::ShapeFunction;
-                _builder[std::type_index(typeid(MeshElement))] =
-                    makeLocalAssemblerBuilder<ShapeFunction>();
-            });
-    }
-
-    /// Returns data pointer to the newly created local assembler data.
-    ///
-    /// \attention
-    /// The index \c id is not necessarily the mesh item's id. Especially when
-    /// having multiple meshes it will differ from the latter.
-    LADataIntfPtr operator()(std::size_t const id,
-                             MeshLib::Element const& mesh_item,
-                             ConstructorArgs&&... args) const
-    {
-        auto const type_idx = std::type_index(typeid(mesh_item));
-        auto const it = _builder.find(type_idx);
-
-        if (it != _builder.end())
-        {
-            auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
-            return it->second(mesh_item, num_local_dof,
-                              std::forward<ConstructorArgs>(args)...);
-        }
-        OGS_FATAL(
-            "You are trying to build a local assembler for an unknown mesh "
-            "element type ({:s})."
-            " Maybe you have disabled this mesh element type in your build "
-            "configuration, or a mesh element order does not match shape "
-            "function order given in the project file.",
-            type_idx.name());
-    }
-
-private:
-    using LADataBuilder =
-        std::function<LADataIntfPtr(MeshLib::Element const& e,
-                                    std::size_t const local_matrix_size,
-                                    ConstructorArgs&&...)>;
-
-    template <typename ShapeFunction>
-    using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
-        typename ShapeFunction::MeshElement>::IntegrationMethod;
-
-    template <typename ShapeFunction>
-    using LAData =
-        LocalAssemblerData<ShapeFunction, IntegrationMethod<ShapeFunction>,
-                           GlobalDim>;
-
-    /// Generates a function that creates a new LocalAssembler of type
-    /// LAData<ShapeFunction>. Only functions with shape function's dimension
-    /// less or equal to the global dimension are instantiated, e.g.  following
-    /// combinations of shape functions and global dimensions: (Line2, 1),
-    /// (Line2, 2), (Line2, 3), (Hex20, 3) but not (Hex20, 2) or (Hex20, 1).
-    template <typename ShapeFunction>
-    static LADataBuilder makeLocalAssemblerBuilder()
-    {
-        return [](MeshLib::Element const& e,
-                  std::size_t const local_matrix_size,
-                  ConstructorArgs&&... args)
-        {
-            return LADataIntfPtr{new LAData<ShapeFunction>{
-                e, local_matrix_size, std::forward<ConstructorArgs>(args)...}};
-        };
-    }
-
-    /// Mapping of element types to local assembler constructors.
-    std::unordered_map<std::type_index, LADataBuilder> _builder;
-
-    NumLib::LocalToGlobalIndexMap const& _dof_table;
-};
-
-}  // namespace ProcessLib
diff --git a/ProcessLib/Utils/LocalDataInitializerHM.h b/ProcessLib/Utils/LocalDataInitializerHM.h
deleted file mode 100644
index ba4eb33b8f4352380470121f5a24623bbe339bcf..0000000000000000000000000000000000000000
--- a/ProcessLib/Utils/LocalDataInitializerHM.h
+++ /dev/null
@@ -1,151 +0,0 @@
-/**
- * \file
- * \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
-
-#include <functional>
-#include <memory>
-#include <type_traits>
-#include <typeindex>
-#include <typeinfo>
-#include <unordered_map>
-
-#include "EnabledElements.h"
-#include "MeshLib/Elements/Elements.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "NumLib/Fem/FiniteElement/LowerDimShapeTable.h"
-#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
-
-namespace ProcessLib
-{
-/// The LocalDataInitializerHM is a functor creating a local assembler data with
-/// corresponding to the mesh element type shape functions and calling
-/// initialization of the new local assembler data.
-/// For example for MeshLib::Quad a local assembler data with template argument
-/// NumLib::ShapeQuad4 is created.
-///
-/// \attention This is modified version of the ProcessLib::LocalDataInitializer
-/// class which does not include line or point elements. For the shape functions
-/// of order 2 (used for displacement) a shape function of order 1 will be used
-/// for the scalar variables.
-template <typename LocalAssemblerInterface,
-          template <typename, typename, typename, int> class LocalAssemblerData,
-          int GlobalDim, typename... ConstructorArgs>
-class LocalDataInitializerHM final
-{
-    struct IsElementEnabled
-    {
-        template <typename ElementTraits>
-        constexpr bool operator()(ElementTraits*) const
-        {
-            if constexpr (GlobalDim < ElementTraits::ShapeFunction::DIM)
-            {
-                return false;
-            }
-
-            // mechanics processes in OGS are defined in 2D and 3D only
-            return ElementTraits::Element::dimension >= 2;
-        }
-    };
-
-public:
-    using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
-
-    explicit LocalDataInitializerHM(
-        NumLib::LocalToGlobalIndexMap const& dof_table)
-        : _dof_table(dof_table)
-    {
-        using EnabledElementTraits =
-            decltype(BaseLib::TMP::filter<EnabledElementTraitsLagrange>(
-                std::declval<IsElementEnabled>()));
-
-        BaseLib::TMP::foreach<EnabledElementTraits>(
-            [this]<typename ET>(ET*)
-            {
-                using MeshElement = typename ET::Element;
-                using Shp = typename ET::ShapeFunction;
-                using ShpLow = typename ET::LowerOrderShapeFunction;
-
-                _builder[std::type_index(typeid(MeshElement))] =
-                    makeLocalAssemblerBuilder<Shp, ShpLow>();
-            });
-    }
-
-    /// Returns data pointer to the newly created local assembler data.
-    ///
-    /// \attention
-    /// The index \c id is not necessarily the mesh item's id. Especially when
-    /// having multiple meshes it will differ from the latter.
-    LADataIntfPtr operator()(std::size_t const id,
-                             MeshLib::Element const& mesh_item,
-                             ConstructorArgs&&... args) const
-    {
-        auto const type_idx = std::type_index(typeid(mesh_item));
-        auto const it = _builder.find(type_idx);
-
-        if (it != _builder.end())
-        {
-            auto const num_local_dof = _dof_table.getNumberOfElementDOF(id);
-            return it->second(mesh_item, num_local_dof,
-                              std::forward<ConstructorArgs>(args)...);
-        }
-        OGS_FATAL(
-            "You are trying to build a local assembler for an unknown mesh "
-            "element type ({:s})."
-            " Maybe you have disabled this mesh element type in your build "
-            "configuration, or a mesh element order does not match shape "
-            "function order given in the project file.",
-            type_idx.name());
-    }
-
-private:
-    using LADataBuilder =
-        std::function<LADataIntfPtr(MeshLib::Element const& e,
-                                    std::size_t const local_matrix_size,
-                                    ConstructorArgs&&...)>;
-
-    template <typename ShapeFunction>
-    using IntegrationMethod = typename NumLib::GaussLegendreIntegrationPolicy<
-        typename ShapeFunction::MeshElement>::IntegrationMethod;
-
-    template <typename ShapeFunctionDisplacement,
-              typename ShapeFunctionPressure>
-    using LAData =
-        LocalAssemblerData<ShapeFunctionDisplacement, ShapeFunctionPressure,
-                           IntegrationMethod<ShapeFunctionDisplacement>,
-                           GlobalDim>;
-
-    /// Generates a function that creates a new LocalAssembler of type
-    /// LAData<ShapeFunctionDisplacement, ShapeFunctionPressure>. Only functions
-    /// with shape function's dimension less or equal to the global dimension
-    /// are instantiated, e.g.  following combinations of shape functions and
-    /// global dimensions: (Line2, 1), (Line2, 2), (Line2, 3), (Hex20, 3) but
-    /// not (Hex20, 2) or (Hex20, 1).
-    template <typename ShapeFunction, typename LowerOrderShapeFunction>
-    static LADataBuilder makeLocalAssemblerBuilder()
-    {
-        return [](MeshLib::Element const& e,
-                  std::size_t const local_matrix_size,
-                  ConstructorArgs&&... args)
-        {
-            return LADataIntfPtr{
-                new LAData<ShapeFunction, LowerOrderShapeFunction>{
-                    e, local_matrix_size,
-                    std::forward<ConstructorArgs>(args)...}};
-        };
-    }
-
-    /// Mapping of element types to local assembler constructors.
-    std::unordered_map<std::type_index, LADataBuilder> _builder;
-
-    NumLib::LocalToGlobalIndexMap const& _dof_table;
-};
-
-}  // namespace ProcessLib
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index 9a9b2dae8ece8028dd0dd9260b85c7f20bd17ef5..aff128d2b652f787da8fdf44ff8713f2c566b918 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -24,7 +24,6 @@
 #include "NumLib/Function/Interpolation.h"
 #include "NumLib/NumericsConfig.h"
 #include "ProcessLib/Utils/CreateLocalAssemblers.h"
-#include "ProcessLib/Utils/LocalDataInitializer.h"
 #include "Tests/VectorUtils.h"
 
 namespace ExtrapolationTest