diff --git a/GeoLib/Grid.h b/GeoLib/Grid.h
index d119ac3e87911eb7645ff097eb086674044569af..46a521ae0c9bc7e918e4f0e814d1a150a606a935 100644
--- a/GeoLib/Grid.h
+++ b/GeoLib/Grid.h
@@ -50,7 +50,7 @@ public:
      */
     template <typename InputIterator>
     Grid(InputIterator first, InputIterator last,
-         std::size_t max_num_per_grid_cell = 512);
+         std::size_t max_num_per_grid_cell = 16);
 
     Grid(Grid const&) = delete;
     Grid& operator=(Grid const&) = delete;
diff --git a/MeshGeoToolsLib/IdentifySubdomainMesh.cpp b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
index 9d210eb60e6a49ae8c851f2e4a0b87a503c06655..99d5012998fcbee9c21675c991be4cef07ae80b0 100644
--- a/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
+++ b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
@@ -7,9 +7,11 @@
  *              http://www.opengeosys.org/project/license
  */
 
-#include <map>
+#include <range/v3/range/conversion.hpp>
+#include <unordered_map>
 #include <vector>
 
+#include "BaseLib/RunTime.h"
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/Node.h"
@@ -48,39 +50,26 @@ std::vector<std::size_t> identifySubdomainMeshNodes(
 /// \note It is recommended to include only base nodes of elements into the
 /// search \c node_ids.
 std::vector<std::size_t> findElementsInMesh(
-    MeshLib::Mesh const& mesh, std::vector<std::size_t> const& node_ids)
+    std::vector<std::size_t> const& node_ids,
+    std::vector<std::vector<std::size_t>> const& connected_element_ids_per_node)
 {
-    //
-    // Collect all element ids for all nodes.
-    //
-    std::vector<std::size_t> common_element_ids;
-    // Every node is connected to at least one element.
-    auto const nnodes = node_ids.size();
-    common_element_ids.reserve(nnodes);
-
-    for (auto const node_id : node_ids)
-    {
-        auto const& connected_elements =
-            mesh.getElementsConnectedToNode(node_id);
-        std::transform(begin(connected_elements), end(connected_elements),
-                       back_inserter(common_element_ids),
-                       [](MeshLib::Element const* const e)
-                       { return e->getID(); });
-    }
-
     //
     // Count how often an element is shared by all nodes.
     //
-    std::map<std::size_t, int> element_counts;
-    for (auto const element_id : common_element_ids)
+    std::unordered_map<std::size_t, int> element_counts(8);
+    for (auto const node_id : node_ids)
     {
-        element_counts[element_id]++;
+        for (auto const element_id : connected_element_ids_per_node[node_id])
+        {
+            element_counts[element_id]++;
+        }
     }
 
     //
     // Elements which are shared by as many nodes as the input nodes are the
     // desired elements.
     //
+    auto const nnodes = node_ids.size();
     std::vector<std::size_t> element_ids;
     for (auto const& pair : element_counts)
     {
@@ -105,6 +94,15 @@ std::vector<std::vector<std::size_t>> identifySubdomainMeshElements(
     std::vector<std::vector<std::size_t>> bulk_element_ids_map(
         subdomain_mesh.getNumberOfElements());
 
+    // For each node a vector of connected element ids of that node.
+    std::vector<std::vector<std::size_t>> connected_element_ids_per_node(
+        bulk_mesh.getNumberOfNodes());
+    for (auto const node_id : bulk_mesh.getNodes() | MeshLib::views::ids)
+    {
+        connected_element_ids_per_node[node_id] =
+            bulk_mesh.getElementsConnectedToNode(node_id) |
+            MeshLib::views::ids | ranges::to<std::vector>;
+    }
     for (auto* const e : subdomain_mesh.getElements())
     {
         std::vector<std::size_t> element_node_ids(e->getNumberOfBaseNodes());
@@ -118,8 +116,9 @@ std::vector<std::vector<std::size_t>> identifySubdomainMeshElements(
                        begin(element_node_ids_bulk),
                        [&bulk_node_ids](std::size_t const id)
                        { return bulk_node_ids[id]; });
-        std::vector<std::size_t> bulk_element_ids =
-            findElementsInMesh(bulk_mesh, element_node_ids_bulk);
+
+        std::vector<std::size_t> bulk_element_ids = findElementsInMesh(
+            element_node_ids_bulk, connected_element_ids_per_node);
 
         if (bulk_element_ids.empty())
         {
@@ -195,15 +194,22 @@ void identifySubdomainMesh(MeshLib::Mesh& subdomain_mesh,
                            MeshNodeSearcher const& mesh_node_searcher,
                            bool const force_overwrite = false)
 {
+    BaseLib::RunTime time;
+    time.start();
     auto const& bulk_node_ids =
         identifySubdomainMeshNodes(subdomain_mesh, mesh_node_searcher);
+    INFO("identifySubdomainMesh(): identifySubdomainMeshNodes took {:g} s",
+         time.elapsed());
 
     updateOrCheckExistingSubdomainProperty(
         subdomain_mesh, MeshLib::getBulkIDString(MeshLib::MeshItemType::Node),
         bulk_node_ids, MeshLib::MeshItemType::Node, force_overwrite);
 
+    time.start();
     auto const& bulk_element_ids =
         identifySubdomainMeshElements(subdomain_mesh, bulk_mesh);
+    INFO("identifySubdomainMesh(): identifySubdomainMeshElements took {:g} s",
+         time.elapsed());
 
     // The bulk_element_ids could be of two types: one element per entry---this
     // is the expected case for the boundary meshes; multiple elements per
diff --git a/MeshLib/Elements/Element.h b/MeshLib/Elements/Element.h
index 30a2b233bd8cb851d63b5645504cc4091f11f56a..400dba1b6108970a1f31625b015ec492c48102c5 100644
--- a/MeshLib/Elements/Element.h
+++ b/MeshLib/Elements/Element.h
@@ -86,7 +86,7 @@ public:
     virtual const Element* getBoundary(unsigned i) const = 0;
 
     /// Returns the ID of the element.
-    virtual std::size_t getID() const final { return _id; }
+    std::size_t getID() const { return _id; }
 
     virtual unsigned getNumberOfBoundaries() const = 0;
 
@@ -195,7 +195,7 @@ protected:
     explicit Element(std::size_t id);
 
     /// Sets the element ID.
-    virtual void setID(std::size_t id) final { _id = id; }
+    void setID(std::size_t id) { _id = id; }
 
     std::size_t _id;