diff --git a/NumLib/DOF/LocalToGlobalIndexMap.cpp b/NumLib/DOF/LocalToGlobalIndexMap.cpp
index c435af4b98211ca7dd2c884bc85c8c30efe0e58a..d3e0166654c4cebcf11455864a37f1b706f2ecc1 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.cpp
+++ b/NumLib/DOF/LocalToGlobalIndexMap.cpp
@@ -98,15 +98,14 @@ void LocalToGlobalIndexMap::findGlobalIndices(
 }
 
 LocalToGlobalIndexMap::LocalToGlobalIndexMap(
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+    std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
     NumLib::ComponentOrder const order)
     : LocalToGlobalIndexMap(std::move(mesh_subsets), std::vector<unsigned>(mesh_subsets.size(), 1), order)
 {
 }
 
-
 LocalToGlobalIndexMap::LocalToGlobalIndexMap(
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+    std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
     std::vector<unsigned> const& vec_var_n_components,
     NumLib::ComponentOrder const order)
     : _mesh_subsets(std::move(mesh_subsets)),
@@ -127,7 +126,7 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
             auto const global_component_id =
                 getGlobalComponent(variable_id, component_id);
 
-            auto const& mss = *_mesh_subsets[global_component_id];
+            auto const& mss = _mesh_subsets[global_component_id];
             for (int subset_id = 0; subset_id < static_cast<int>(mss.size());
                  ++subset_id)
             {
@@ -143,11 +142,10 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
     }
 }
 
-
 LocalToGlobalIndexMap::LocalToGlobalIndexMap(
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+    std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
     std::vector<unsigned> const& vec_var_n_components,
-    std::vector<std::vector<MeshLib::Element*>const*> const& vec_var_elements,
+    std::vector<std::vector<MeshLib::Element*> const*> const& vec_var_elements,
     NumLib::ComponentOrder const order)
     : _mesh_subsets(std::move(mesh_subsets)),
       _mesh_component_map(_mesh_subsets, order),
@@ -178,7 +176,7 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
             auto const global_component_id =
                 getGlobalComponent(variable_id, component_id);
 
-            auto const& mss = *_mesh_subsets[global_component_id];
+            auto const& mss = _mesh_subsets[global_component_id];
             for (int subset_id = 0; subset_id < static_cast<int>(mss.size());
                  ++subset_id)
             {
@@ -195,7 +193,7 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
 }
 
 LocalToGlobalIndexMap::LocalToGlobalIndexMap(
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+    std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
     std::vector<int> const& global_component_ids,
     std::vector<MeshLib::Element*> const& elements,
     NumLib::MeshComponentMap&& mesh_component_map)
@@ -213,7 +211,7 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
 
     for (int i = 0; i < global_component_ids.size(); ++i)
     {
-        auto const& mss = *_mesh_subsets[i];
+        auto const& mss = _mesh_subsets[i];
 
         // For all MeshSubset in mesh_subsets and each element of that
         // MeshSubset
@@ -232,7 +230,7 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
 LocalToGlobalIndexMap* LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
     int const variable_id,
     std::vector<int> const& component_ids,
-    std::unique_ptr<MeshLib::MeshSubsets>&& mesh_subsets,
+    MeshLib::MeshSubsets&& mesh_subsets,
     std::vector<MeshLib::Element*> const& elements) const
 {
     DBUG("Construct reduced local to global index map.");
@@ -248,13 +246,13 @@ LocalToGlobalIndexMap* LocalToGlobalIndexMap::deriveBoundaryConstrainedMap(
             getGlobalComponent(variable_id, component_id));
 
     auto mesh_component_map =
-        _mesh_component_map.getSubset(global_component_ids, *mesh_subsets);
+        _mesh_component_map.getSubset(global_component_ids, mesh_subsets);
 
     // Create copies of the mesh_subsets for each of the global components.
     // The last component is moved after the for-loop.
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
+    std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
     for (int i = 0; i < static_cast<int>(global_component_ids.size()) - 1; ++i)
-        all_mesh_subsets.emplace_back(new MeshLib::MeshSubsets{*mesh_subsets});
+        all_mesh_subsets.emplace_back(mesh_subsets);
     all_mesh_subsets.emplace_back(std::move(mesh_subsets));
 
     return new LocalToGlobalIndexMap(std::move(all_mesh_subsets),
diff --git a/NumLib/DOF/LocalToGlobalIndexMap.h b/NumLib/DOF/LocalToGlobalIndexMap.h
index 8bfadc2191e1d53d30b820c725c660ca71cb7727..c14202927978003b62e0a2207d509b1bb5063a04 100644
--- a/NumLib/DOF/LocalToGlobalIndexMap.h
+++ b/NumLib/DOF/LocalToGlobalIndexMap.h
@@ -47,7 +47,7 @@ public:
     /// \attention This constructor assumes the number of the given mesh subsets
     /// is equal to the number of variables, i.e. every variable has a single component.
     LocalToGlobalIndexMap(
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+        std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
         NumLib::ComponentOrder const order);
 
     /// Creates a MeshComponentMap internally and stores the global indices for
@@ -59,7 +59,7 @@ public:
     /// should be equal to the size of the mesh_subsets.
     /// \param order  type of ordering values in a vector
     LocalToGlobalIndexMap(
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+        std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
         std::vector<unsigned> const& vec_var_n_components,
         NumLib::ComponentOrder const order);
 
@@ -73,7 +73,7 @@ public:
     /// \param vec_var_elements  a vector of active mesh elements for each variable.
     /// \param order  type of ordering values in a vector
     LocalToGlobalIndexMap(
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+        std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
         std::vector<unsigned> const& vec_var_n_components,
         std::vector<std::vector<MeshLib::Element*>const*> const& vec_var_elements,
         NumLib::ComponentOrder const order);
@@ -86,7 +86,7 @@ public:
     LocalToGlobalIndexMap* deriveBoundaryConstrainedMap(
         int const variable_id,
         std::vector<int> const& component_ids,
-        std::unique_ptr<MeshLib::MeshSubsets>&& mesh_subsets,
+        MeshLib::MeshSubsets&& mesh_subsets,
         std::vector<MeshLib::Element*> const& elements) const;
 
     /// Returns total number of degrees of freedom including those located in
@@ -159,7 +159,7 @@ public:
 
     MeshLib::MeshSubsets const& getMeshSubsets(int const global_component_id) const
     {
-        return *_mesh_subsets[global_component_id];
+        return _mesh_subsets[global_component_id];
     }
 
 private:
@@ -169,7 +169,7 @@ private:
     /// \attention The passed mesh_component_map is in undefined state after
     /// this construtor.
     explicit LocalToGlobalIndexMap(
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
+        std::vector<MeshLib::MeshSubsets>&& mesh_subsets,
         std::vector<int> const& global_component_ids,
         std::vector<MeshLib::Element*> const& elements,
         NumLib::MeshComponentMap&& mesh_component_map);
@@ -197,7 +197,7 @@ private:
 
 private:
     /// A vector of mesh subsets for each process variables' components.
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> const _mesh_subsets;
+    std::vector<MeshLib::MeshSubsets> const _mesh_subsets;
     NumLib::MeshComponentMap _mesh_component_map;
 
     using Table = Eigen::Matrix<LineIndex, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
diff --git a/NumLib/DOF/MeshComponentMap.cpp b/NumLib/DOF/MeshComponentMap.cpp
index 700767359373bd8e48c5448dd9d90ea530e235b8..19a4b01c701b95e2c59959a69346c5908c9d2a11 100644
--- a/NumLib/DOF/MeshComponentMap.cpp
+++ b/NumLib/DOF/MeshComponentMap.cpp
@@ -28,19 +28,17 @@ GlobalIndexType const MeshComponentMap::nop =
 
 #ifdef USE_PETSC
 MeshComponentMap::MeshComponentMap(
-    const std::vector<std::unique_ptr<MeshLib::MeshSubsets>>& components,
-    ComponentOrder order)
+    const std::vector<MeshLib::MeshSubsets>& components, ComponentOrder order)
 {
     // get number of unknows
     GlobalIndexType num_unknowns = 0;
     for (auto const& c : components)
     {
-        assert(c != nullptr);
-        for (unsigned mesh_subset_index = 0; mesh_subset_index < c->size();
+        for (unsigned mesh_subset_index = 0; mesh_subset_index < c.size();
              mesh_subset_index++)
         {
             MeshLib::MeshSubset const& mesh_subset =
-                c->getMeshSubset(mesh_subset_index);
+                c.getMeshSubset(mesh_subset_index);
             // PETSc always works with MeshLib::NodePartitionedMesh.
             const MeshLib::NodePartitionedMesh& mesh =
                 static_cast<const MeshLib::NodePartitionedMesh&>(
@@ -56,12 +54,11 @@ MeshComponentMap::MeshComponentMap(
     _num_local_dof = 0;
     for (auto const& c : components)
     {
-        assert(c != nullptr);
-        for (unsigned mesh_subset_index = 0; mesh_subset_index < c->size();
+        for (unsigned mesh_subset_index = 0; mesh_subset_index < c.size();
              mesh_subset_index++)
         {
             MeshLib::MeshSubset const& mesh_subset =
-                c->getMeshSubset(mesh_subset_index);
+                c.getMeshSubset(mesh_subset_index);
             assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(
                        &mesh_subset.getMesh()) != nullptr);
             std::size_t const mesh_id = mesh_subset.getMeshID();
@@ -117,18 +114,16 @@ MeshComponentMap::MeshComponentMap(
 }
 #else
 MeshComponentMap::MeshComponentMap(
-    const std::vector<std::unique_ptr<MeshLib::MeshSubsets>>& components,
-    ComponentOrder order)
+    const std::vector<MeshLib::MeshSubsets>& components, ComponentOrder order)
 {
     // construct dict (and here we number global_index by component type)
     GlobalIndexType global_index = 0;
     std::size_t comp_id = 0;
     for (auto const& c : components)
     {
-        assert (c != nullptr);
-        for (std::size_t mesh_subset_index = 0; mesh_subset_index < c->size(); mesh_subset_index++)
+        for (std::size_t mesh_subset_index = 0; mesh_subset_index < c.size(); mesh_subset_index++)
         {
-            MeshLib::MeshSubset const& mesh_subset = c->getMeshSubset(mesh_subset_index);
+            MeshLib::MeshSubset const& mesh_subset = c.getMeshSubset(mesh_subset_index);
             std::size_t const mesh_id = mesh_subset.getMeshID();
             // mesh items are ordered first by node, cell, ....
             for (std::size_t j=0; j<mesh_subset.getNumberOfNodes(); j++)
diff --git a/NumLib/DOF/MeshComponentMap.h b/NumLib/DOF/MeshComponentMap.h
index e4402cc7813055b34479cb47822852259cbed6cd..7bc4c9618d828deb53f0461d3693aa9c05b1850b 100644
--- a/NumLib/DOF/MeshComponentMap.h
+++ b/NumLib/DOF/MeshComponentMap.h
@@ -38,9 +38,8 @@ public:
 public:
     /// \param components   a vector of components
     /// \param order        type of ordering values in a vector
-    MeshComponentMap(
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>> const& components,
-        ComponentOrder order);
+    MeshComponentMap(std::vector<MeshLib::MeshSubsets> const& components,
+                     ComponentOrder order);
 
     /// Creates a multi-component subset of the current mesh component map.
     /// The order (BY_LOCATION/BY_COMPONENT) of components is the same as of the
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index fe4949de469d88ac8051f343ac94978f5bfec567..d2f0ea3ea4c5086546f44ce24c59bc1b5701324b 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -49,8 +49,7 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
     // to each of the MeshSubset in the mesh_subsets.
     _mesh_subset_all_nodes.reset(
         mesh_subsets.getMeshSubset(0).getIntersectionByNodes(nodes));
-    std::unique_ptr<MeshLib::MeshSubsets> all_mesh_subsets{
-        new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}};
+    MeshLib::MeshSubsets all_mesh_subsets{_mesh_subset_all_nodes.get()};
 
     // Create local DOF table from intersected mesh subsets for the given
     // variable and component ids.
@@ -58,8 +57,8 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
         variable_id, {component_id}, std::move(all_mesh_subsets), _elements));
 
     createLocalAssemblers<LocalAssemblerImplementation>(
-        global_dim, _elements, *_dof_table_boundary, shapefunction_order, _local_assemblers,
-        is_axially_symmetric, _integration_order, _data);
+        global_dim, _elements, *_dof_table_boundary, shapefunction_order,
+        _local_assemblers, is_axially_symmetric, _integration_order, _data);
 }
 
 template <typename BoundaryConditionData,
diff --git a/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp b/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp
index c6f31f5e5505253afa384b86c0047640ad7c8242..6dba93cf506f693c296f2ce11f4ee3d969b51477 100644
--- a/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp
+++ b/ProcessLib/CalculateSurfaceFlux/CalculateSurfaceFlux.cpp
@@ -29,14 +29,11 @@ CalculateSurfaceFlux::CalculateSurfaceFlux(
         new MeshLib::MeshSubset(boundary_mesh, &boundary_mesh.getNodes()));
 
     // Collect the mesh subsets in a vector.
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
+    std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
     std::generate_n(
         std::back_inserter(all_mesh_subsets),
         bulk_property_number_of_components,
-        [&]() {
-            return std::unique_ptr<MeshLib::MeshSubsets>{
-                new MeshLib::MeshSubsets{mesh_subset_all_nodes.get()}};
-        });
+        [&]() { return MeshLib::MeshSubsets{mesh_subset_all_nodes.get()}; });
 
     // needed for creation of local assemblers
     std::unique_ptr<NumLib::LocalToGlobalIndexMap const> dof_table(
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
index 4496c7fa41c5e4ef114593806b0d74cab0ca0eb2..4ca7953b5d4988dbe65d94def06bb698e6e0a599 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsProcess.h
@@ -71,17 +71,15 @@ private:
         // Collect the mesh subsets in a vector.
 
         // For pressure, which is the first
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
-        all_mesh_subsets.push_back(std::unique_ptr<MeshLib::MeshSubsets>{
-            new MeshLib::MeshSubsets{_mesh_subset_base_nodes.get()}});
+        std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
+        all_mesh_subsets.emplace_back(_mesh_subset_base_nodes.get());
 
         // For displacement.
         std::generate_n(
             std::back_inserter(all_mesh_subsets),
             getProcessVariables()[1].get().getNumberOfComponents(),
             [&]() {
-                return std::unique_ptr<MeshLib::MeshSubsets>{
-                    new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}};
+                return MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()};
             });
 
         std::vector<unsigned> const vec_n_components{1, DisplacementDim};
@@ -105,10 +103,9 @@ private:
 
         // TODO move the two data members somewhere else.
         // for extrapolation of secondary variables
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>>
-            all_mesh_subsets_single_component;
+        std::vector<MeshLib::MeshSubsets> all_mesh_subsets_single_component;
         all_mesh_subsets_single_component.emplace_back(
-            new MeshLib::MeshSubsets(_mesh_subset_all_nodes.get()));
+            _mesh_subset_all_nodes.get());
         _local_to_global_index_map_single_component.reset(
             new NumLib::LocalToGlobalIndexMap(
                 std::move(all_mesh_subsets_single_component),
diff --git a/ProcessLib/LIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h b/ProcessLib/LIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
index 0509fe0d2427ebb6686d725af2c1da649f0bada5..54a41e4a48f79ad64f6aec7a15baeba2eb31124c 100644
--- a/ProcessLib/LIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
+++ b/ProcessLib/LIE/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
@@ -56,8 +56,7 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
     // to each of the MeshSubset in the mesh_subsets.
     _mesh_subset_all_nodes.reset(
         mesh_subsets.getMeshSubset(0).getIntersectionByNodes(nodes));
-    std::unique_ptr<MeshLib::MeshSubsets> all_mesh_subsets{
-        new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}};
+    MeshLib::MeshSubsets all_mesh_subsets{_mesh_subset_all_nodes.get()};
 
     // Create local DOF table from intersected mesh subsets for the given
     // variable and component ids.
@@ -65,8 +64,9 @@ GenericNaturalBoundaryCondition<BoundaryConditionData,
         variable_id, {component_id}, std::move(all_mesh_subsets), _elements));
 
     createLocalAssemblers<LocalAssemblerImplementation>(
-        global_dim, _elements, *_dof_table_boundary, shapefunction_order, _local_assemblers,
-        is_axially_symmetric, _integration_order, _data, fracture_prop, variable_id);
+        global_dim, _elements, *_dof_table_boundary, shapefunction_order,
+        _local_assemblers, is_axially_symmetric, _integration_order, _data,
+        fracture_prop, variable_id);
 }
 
 template <typename BoundaryConditionData,
diff --git a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
index 9b6be35e0d9f86a0d55cf348eb04422d986968d0..d9ba570241f02c83eabc2bb3925c3562d36972d3 100644
--- a/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
+++ b/ProcessLib/LIE/HydroMechanics/HydroMechanicsProcess.cpp
@@ -143,13 +143,12 @@ void HydroMechanicsProcess<GlobalDim>::constructDofTable()
     }
 
     // Collect the mesh subsets in a vector.
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
+    std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
     std::vector<unsigned> vec_n_components;
     std::vector<std::vector<MeshLib::Element*> const*> vec_var_elements;
     // pressure
     vec_n_components.push_back(1);
-    all_mesh_subsets.push_back(std::unique_ptr<MeshLib::MeshSubsets>{
-        new MeshLib::MeshSubsets{_mesh_subset_nodes_p.get()}});
+    all_mesh_subsets.emplace_back(_mesh_subset_nodes_p.get());
     if (!_process_data.deactivate_matrix_in_flow)
     {
         vec_var_elements.push_back(&_mesh.getElements());
@@ -163,8 +162,7 @@ void HydroMechanicsProcess<GlobalDim>::constructDofTable()
     // regular displacement
     vec_n_components.push_back(GlobalDim);
     std::generate_n(std::back_inserter(all_mesh_subsets), GlobalDim, [&]() {
-        return std::unique_ptr<MeshLib::MeshSubsets>{
-            new MeshLib::MeshSubsets{_mesh_subset_matrix_nodes.get()}};
+        return MeshLib::MeshSubsets{_mesh_subset_matrix_nodes.get()};
     });
     vec_var_elements.push_back(&_vec_matrix_elements);
     if (!_vec_fracture_nodes.empty())
@@ -172,8 +170,7 @@ void HydroMechanicsProcess<GlobalDim>::constructDofTable()
         // displacement jump
         vec_n_components.push_back(GlobalDim);
         std::generate_n(std::back_inserter(all_mesh_subsets), GlobalDim, [&]() {
-            return std::unique_ptr<MeshLib::MeshSubsets>{
-                new MeshLib::MeshSubsets{_mesh_subset_fracture_nodes.get()}};
+            return MeshLib::MeshSubsets{_mesh_subset_fracture_nodes.get()};
         });
         vec_var_elements.push_back(&_vec_fracture_matrix_elements);
     }
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
index ecf7bc3fc81521e3b235101c8cb11598bfa92e82..9be631e3a50c48c589996cf36a73b1af73f33637 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
@@ -125,28 +125,20 @@ void SmallDeformationProcess<DisplacementDim>::constructDofTable()
     {
         _mesh_subset_fracture_nodes.push_back(
             std::unique_ptr<MeshLib::MeshSubset const>(
-                        new MeshLib::MeshSubset(_mesh, &_vec_fracture_nodes[i])
-                        ));
+                new MeshLib::MeshSubset(_mesh, &_vec_fracture_nodes[i])));
     }
 
     // Collect the mesh subsets in a vector.
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
+    std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
     std::generate_n(
-        std::back_inserter(all_mesh_subsets),
-        DisplacementDim,
-        [&]() {
-            return std::unique_ptr<MeshLib::MeshSubsets>{
-                new MeshLib::MeshSubsets{_mesh_subset_matrix_nodes.get()}};
+        std::back_inserter(all_mesh_subsets), DisplacementDim, [&]() {
+            return MeshLib::MeshSubsets{_mesh_subset_matrix_nodes.get()};
         });
     for (auto& ms : _mesh_subset_fracture_nodes)
     {
-        std::generate_n(
-            std::back_inserter(all_mesh_subsets),
-            DisplacementDim,
-            [&]() {
-                return std::unique_ptr<MeshLib::MeshSubsets>{
-                    new MeshLib::MeshSubsets{ms.get()}};
-            });
+        std::generate_n(std::back_inserter(all_mesh_subsets),
+                        DisplacementDim,
+                        [&]() { return MeshLib::MeshSubsets{ms.get()}; });
     }
 
     std::vector<unsigned> const vec_n_components(
@@ -183,10 +175,9 @@ void SmallDeformationProcess<DisplacementDim>::initializeConcreteProcess(
 
     // TODO move the two data members somewhere else.
     // for extrapolation of secondary variables
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>>
-        all_mesh_subsets_single_component;
+    std::vector<MeshLib::MeshSubsets> all_mesh_subsets_single_component;
     all_mesh_subsets_single_component.emplace_back(
-        new MeshLib::MeshSubsets(_mesh_subset_all_nodes.get()));
+        _mesh_subset_all_nodes.get());
     _local_to_global_index_map_single_component.reset(
         new NumLib::LocalToGlobalIndexMap(
             std::move(all_mesh_subsets_single_component),
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 0fb1ca219366f8b317323026742b4c97d6404042..4fc07febe18116d4d412b4fb1af0479a8d1012e0 100644
--- a/ProcessLib/Process.cpp
+++ b/ProcessLib/Process.cpp
@@ -148,15 +148,14 @@ void Process::constructDofTable()
         new MeshLib::MeshSubset(_mesh, &_mesh.getNodes()));
 
     // Collect the mesh subsets in a vector.
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
+    std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
     for (ProcessVariable const& pv : _process_variables)
     {
         std::generate_n(
             std::back_inserter(all_mesh_subsets),
             pv.getNumberOfComponents(),
             [&]() {
-                return std::unique_ptr<MeshLib::MeshSubsets>{
-                    new MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()}};
+                return MeshLib::MeshSubsets{_mesh_subset_all_nodes.get()};
             });
     }
 
@@ -186,10 +185,9 @@ void Process::initializeExtrapolator()
     else
     {
         // Otherwise construct a new DOF table.
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>>
-            all_mesh_subsets_single_component;
+        std::vector<MeshLib::MeshSubsets> all_mesh_subsets_single_component;
         all_mesh_subsets_single_component.emplace_back(
-            new MeshLib::MeshSubsets(_mesh_subset_all_nodes.get()));
+            _mesh_subset_all_nodes.get());
 
         dof_table_single_component = new NumLib::LocalToGlobalIndexMap(
             std::move(all_mesh_subsets_single_component),
diff --git a/ProcessLib/SmallDeformation/SmallDeformationProcess.h b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
index ff536b335f3904323749462c92748d29cc88e8c8..1effe277d92f031abe5100c4a5a417a2cb46940f 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationProcess.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationProcess.h
@@ -71,10 +71,9 @@ private:
 
         // TODO move the two data members somewhere else.
         // for extrapolation of secondary variables
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>>
-            all_mesh_subsets_single_component;
+        std::vector<MeshLib::MeshSubsets> all_mesh_subsets_single_component;
         all_mesh_subsets_single_component.emplace_back(
-            new MeshLib::MeshSubsets(_mesh_subset_all_nodes.get()));
+            _mesh_subset_all_nodes.get());
         _local_to_global_index_map_single_component.reset(
             new NumLib::LocalToGlobalIndexMap(
                 std::move(all_mesh_subsets_single_component),
diff --git a/Tests/NumLib/LocalToGlobalIndexMap.cpp b/Tests/NumLib/LocalToGlobalIndexMap.cpp
index 11ce39c3e04fc23bf16c840aef674218d2006324..a4d630081e26e2de857698a5dc5fdab377c20caf 100644
--- a/Tests/NumLib/LocalToGlobalIndexMap.cpp
+++ b/Tests/NumLib/LocalToGlobalIndexMap.cpp
@@ -30,8 +30,8 @@ public:
         nodesSubset.reset(new MeshLib::MeshSubset(*mesh, &mesh->getNodes()));
 
         // Add two components both based on the same nodesSubset.
-        components.emplace_back(new MeshLib::MeshSubsets{nodesSubset.get()});
-        components.emplace_back(new MeshLib::MeshSubsets{nodesSubset.get()});
+        components.emplace_back(nodesSubset.get());
+        components.emplace_back(nodesSubset.get());
     }
 
 protected:
@@ -42,7 +42,7 @@ protected:
     //data component 0 and 1 are assigned to all nodes in the mesh
     static std::size_t const comp0_id = 0;
     static std::size_t const comp1_id = 1;
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> components;
+    std::vector<MeshLib::MeshSubsets> components;
 
     std::unique_ptr<NumLib::LocalToGlobalIndexMap const> dof_map;
 };
@@ -103,8 +103,7 @@ TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_SubsetByComponent)
 
     auto const selected_subset = std::unique_ptr<MeshLib::MeshSubset const>{
         nodesSubset->getIntersectionByNodes(selected_nodes)};
-    auto selected_component = std::unique_ptr<MeshLib::MeshSubsets>{
-        new MeshLib::MeshSubsets{selected_subset.get()}};
+    MeshLib::MeshSubsets selected_component{selected_subset.get()};
 
     auto dof_map_subset = std::unique_ptr<NumLib::LocalToGlobalIndexMap>{
         dof_map->deriveBoundaryConstrainedMap(1,    // variable id
@@ -124,7 +123,7 @@ TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_MultipleVariablesMultipleCompon
 #endif
 {
     // test 2 variables (1st variable with 1 component, 2nd variable with 2 components)
-    components.emplace_back(new MeshLib::MeshSubsets{nodesSubset.get()});
+    components.emplace_back(nodesSubset.get());
 
     std::vector<unsigned> vec_var_n_components{1, 2};
 
@@ -153,7 +152,7 @@ TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_MultipleVariablesMultipleCompon
 #endif
 {
     // test 2 variables (1st variable with 2 component, 2nd variable with 1 components)
-    components.emplace_back(new MeshLib::MeshSubsets{nodesSubset.get()});
+    components.emplace_back(nodesSubset.get());
 
     std::vector<unsigned> vec_var_n_components{2, 1};
 
@@ -187,7 +186,7 @@ TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_MultipleVariablesMultipleCompon
     // - 2nd variable with 1 component for nodes of element id 1
     std::vector<MeshLib::Node*> var2_nodes{const_cast<MeshLib::Node*>(mesh->getNode(1)), const_cast<MeshLib::Node*>(mesh->getNode(2))};
     std::unique_ptr<MeshLib::MeshSubset> var2_subset(new MeshLib::MeshSubset(*mesh, &var2_nodes));
-    components.emplace_back(new MeshLib::MeshSubsets{var2_subset.get()});
+    components.emplace_back(var2_subset.get());
 
     std::vector<unsigned> vec_var_n_components{2, 1};
     std::vector<std::vector<MeshLib::Element*>const*> vec_var_elements;
@@ -239,7 +238,7 @@ TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_MultipleVariablesMultipleCompon
     // - 2nd variable with 1 component for 1st node of element id 1
     std::vector<MeshLib::Node*> var2_nodes{const_cast<MeshLib::Node*>(mesh->getNode(1))};
     std::unique_ptr<MeshLib::MeshSubset> var2_subset(new MeshLib::MeshSubset(*mesh, &var2_nodes));
-    components.emplace_back(new MeshLib::MeshSubsets{var2_subset.get()});
+    components.emplace_back(var2_subset.get());
 
     std::vector<unsigned> vec_var_n_components{2, 1};
     std::vector<std::vector<MeshLib::Element*>const*> vec_var_elements;
diff --git a/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp b/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
index 908d622f288a0fc514868759a43418e16196c8b4..df08bff74c27b07decc13e4df44a92dca9cacf71 100644
--- a/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
+++ b/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
@@ -83,18 +83,16 @@ public:
     {
         assert(selected_component < num_components);
 
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>> components;
+        std::vector<MeshLib::MeshSubsets> components;
         for (unsigned i=0; i<num_components; ++i)
         {
-            components.emplace_back(
-                new MeL::MeshSubsets{mesh_items_all_nodes.get()});
+            components.emplace_back(mesh_items_all_nodes.get());
         }
         std::vector<unsigned> vec_var_n_components(1, num_components);
         dof_map.reset(
             new NL::LocalToGlobalIndexMap(std::move(components), vec_var_n_components, order));
 
-        auto components_boundary = std::unique_ptr<MeshLib::MeshSubsets>{
-            new MeL::MeshSubsets{mesh_items_boundary.get()}};
+        MeL::MeshSubsets components_boundary{mesh_items_boundary.get()};
 
         dof_map_boundary.reset(dof_map->deriveBoundaryConstrainedMap(
             0,  // variable id
@@ -110,18 +108,16 @@ public:
     {
         assert(selected_components.size() <= num_components);
 
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>> components;
+        std::vector<MeshLib::MeshSubsets> components;
         for (unsigned i = 0; i < num_components; ++i)
         {
-            components.emplace_back(
-                new MeL::MeshSubsets{mesh_items_all_nodes.get()});
+            components.emplace_back(mesh_items_all_nodes.get());
         }
         std::vector<unsigned> vec_var_n_components(1, num_components);
         dof_map.reset(new NL::LocalToGlobalIndexMap(
             std::move(components), vec_var_n_components, order));
 
-        auto components_boundary = std::unique_ptr<MeshLib::MeshSubsets>{
-            new MeL::MeshSubsets{mesh_items_boundary.get()}};
+        MeL::MeshSubsets components_boundary{mesh_items_boundary.get()};
 
         dof_map_boundary.reset(dof_map->deriveBoundaryConstrainedMap(
             0,  // variable id
diff --git a/Tests/NumLib/TestComponentNorms.cpp b/Tests/NumLib/TestComponentNorms.cpp
index 1d25ec27cad8273c8a7dedc572cc45b957bf05b3..8aff0d9e53cf7de3d410a344884e9268a698f72c 100644
--- a/Tests/NumLib/TestComponentNorms.cpp
+++ b/Tests/NumLib/TestComponentNorms.cpp
@@ -20,13 +20,12 @@
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "Tests/VectorUtils.h"
 
-static std::vector<std::unique_ptr<MeshLib::MeshSubsets>> createMeshSubsets(
+static std::vector<MeshLib::MeshSubsets> createMeshSubsets(
     MeshLib::MeshSubset const& mesh_subset, unsigned const num_components)
 {
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> mesh_subsets;
+    std::vector<MeshLib::MeshSubsets> mesh_subsets;
     for (unsigned i=0; i<num_components; ++i)
-        mesh_subsets.emplace_back(
-            new MeshLib::MeshSubsets{&mesh_subset});
+        mesh_subsets.emplace_back(&mesh_subset);
 
     return mesh_subsets;
 }
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index f9f4ebfc018678ef189c257f9474957189c1668a..a85b1abdb3f2e1a9d94b4757db79ad34835e7fc7 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -135,9 +135,8 @@ public:
         : _integration_order(integration_order)
         , _mesh_subset_all_nodes(mesh, &mesh.getNodes())
     {
-        std::vector<std::unique_ptr<MeshLib::MeshSubsets>> all_mesh_subsets;
-        all_mesh_subsets.emplace_back(
-            new MeshLib::MeshSubsets{&_mesh_subset_all_nodes});
+        std::vector<MeshLib::MeshSubsets> all_mesh_subsets;
+        all_mesh_subsets.emplace_back(&_mesh_subset_all_nodes);
 
         _dof_table.reset(new NumLib::LocalToGlobalIndexMap(
             std::move(all_mesh_subsets), NumLib::ComponentOrder::BY_COMPONENT));
diff --git a/Tests/NumLib/TestMeshComponentMap.cpp b/Tests/NumLib/TestMeshComponentMap.cpp
index 469a8229fb704f10c577e46451df05ca1fa09455..4109be8b5f0d1bf05aca1145af8b2de5edcc713d 100644
--- a/Tests/NumLib/TestMeshComponentMap.cpp
+++ b/Tests/NumLib/TestMeshComponentMap.cpp
@@ -34,8 +34,8 @@ class NumLibMeshComponentMapTest : public ::testing::Test
         nodesSubset = new MeshLib::MeshSubset(*mesh, &mesh->getNodes());
 
         // Add two components both based on the same nodesSubset.
-        components.emplace_back(new MeshLib::MeshSubsets{nodesSubset});
-        components.emplace_back(new MeshLib::MeshSubsets{nodesSubset});
+        components.emplace_back(nodesSubset);
+        components.emplace_back(nodesSubset);
     }
 
     ~NumLibMeshComponentMapTest() override
@@ -52,7 +52,7 @@ class NumLibMeshComponentMapTest : public ::testing::Test
     //data component 0 and 1 are assigned to all nodes in the mesh
     static std::size_t const comp0_id = 0;
     static std::size_t const comp1_id = 1;
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> components;
+    std::vector<MeshLib::MeshSubsets> components;
     MeshComponentMap const* cmap;
 
     //
diff --git a/Tests/NumLib/TestSerialLinearSolver.cpp b/Tests/NumLib/TestSerialLinearSolver.cpp
index a861a895eff8ac53845d6337799cc3159139c006..d929a3f0f8348eb04235736669beb7fcd878d46a 100644
--- a/Tests/NumLib/TestSerialLinearSolver.cpp
+++ b/Tests/NumLib/TestSerialLinearSolver.cpp
@@ -53,8 +53,8 @@ TEST(NumLibSerialLinearSolver, Steady2DdiffusionQuadElem)
     // Allocate a coefficient matrix, RHS and solution vectors
     //--------------------------------------------------------------------------
     // define a mesh item composition in a vector
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> vec_comp_dis;
-    vec_comp_dis.emplace_back(new MeshLib::MeshSubsets{&mesh_items_all_nodes});
+    std::vector<MeshLib::MeshSubsets> vec_comp_dis;
+    vec_comp_dis.emplace_back(&mesh_items_all_nodes);
     NumLib::LocalToGlobalIndexMap local_to_global_index_map(
         std::move(vec_comp_dis), NumLib::ComponentOrder::BY_COMPONENT);
 
diff --git a/Tests/NumLib/TestSparsityPattern.cpp b/Tests/NumLib/TestSparsityPattern.cpp
index 4069824e77f95dfe8c878fa0ca659457379ac067..d03224107d3336d939ff672b3376746a29de5fe5 100644
--- a/Tests/NumLib/TestSparsityPattern.cpp
+++ b/Tests/NumLib/TestSparsityPattern.cpp
@@ -29,8 +29,8 @@ TEST(NumLib_SparsityPattern, DISABLED_SingleComponentLinearMesh)
     std::unique_ptr<MeshLib::MeshSubset const> nodesSubset(
         new MeshLib::MeshSubset(*mesh, &mesh->getNodes()));
 
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> components;
-    components.emplace_back(new MeshLib::MeshSubsets{nodesSubset.get()});
+    std::vector<MeshLib::MeshSubsets> components;
+    components.emplace_back(nodesSubset.get());
     NumLib::LocalToGlobalIndexMap dof_map(
                       std::move(components),
                       NumLib::ComponentOrder::BY_COMPONENT);
@@ -58,8 +58,8 @@ TEST(NumLib_SparsityPattern, DISABLED_SingleComponentQuadraticMesh)
     std::unique_ptr<MeshLib::MeshSubset const> nodesSubset(
         new MeshLib::MeshSubset(*mesh, &mesh->getNodes()));
 
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> components;
-    components.emplace_back(new MeshLib::MeshSubsets{nodesSubset.get()});
+    std::vector<MeshLib::MeshSubsets> components;
+    components.emplace_back(nodesSubset.get());
     NumLib::LocalToGlobalIndexMap dof_map(
                       std::move(components),
                       NumLib::ComponentOrder::BY_COMPONENT);
@@ -88,9 +88,9 @@ TEST(NumLib_SparsityPattern, DISABLED_MultipleComponentsLinearMesh)
     std::unique_ptr<MeshLib::MeshSubset const> nodesSubset(
         new MeshLib::MeshSubset(*mesh, &mesh->getNodes()));
 
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> components;
-    components.emplace_back(new MeshLib::MeshSubsets{nodesSubset.get()});
-    components.emplace_back(new MeshLib::MeshSubsets{nodesSubset.get()});
+    std::vector<MeshLib::MeshSubsets> components;
+    components.emplace_back(nodesSubset.get());
+    components.emplace_back(nodesSubset.get());
     NumLib::LocalToGlobalIndexMap dof_map(
                       std::move(components),
                       NumLib::ComponentOrder::BY_COMPONENT);
@@ -124,9 +124,9 @@ TEST(NumLib_SparsityPattern, DISABLED_MultipleComponentsLinearQuadraticMesh)
     std::unique_ptr<MeshLib::MeshSubset const> allNodesSubset(
         new MeshLib::MeshSubset(*mesh, &mesh->getNodes()));
 
-    std::vector<std::unique_ptr<MeshLib::MeshSubsets>> components;
-    components.emplace_back(new MeshLib::MeshSubsets{baseNodesSubset.get()});
-    components.emplace_back(new MeshLib::MeshSubsets{allNodesSubset.get()});
+    std::vector<MeshLib::MeshSubsets> components;
+    components.emplace_back(baseNodesSubset.get());
+    components.emplace_back(allNodesSubset.get());
     NumLib::LocalToGlobalIndexMap dof_map(
                       std::move(components),
                       NumLib::ComponentOrder::BY_COMPONENT);