From 08dca10c887311b83c80df5f1c68ba8d986cce59 Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Thu, 18 Jan 2024 10:16:28 +0100
Subject: [PATCH] [MGTL] OpenMP parallelization of
 MeshNodeSearcher::getMeshNodeIDs()

It is important to insert the node ids in the same sequence as the points are given.
Therefore, it is necessary to create thread local vectors and insert the content of
the thread local vectors into the node_ids vector with the clause ordered in the
omp directive.
---
 MeshGeoToolsLib/IdentifySubdomainMesh.cpp |  5 ++++-
 MeshGeoToolsLib/MeshNodeSearcher.cpp      | 11 ++++++-----
 2 files changed, 10 insertions(+), 6 deletions(-)

diff --git a/MeshGeoToolsLib/IdentifySubdomainMesh.cpp b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
index 1fdf0c23bca..1ddbdc28ade 100644
--- a/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
+++ b/MeshGeoToolsLib/IdentifySubdomainMesh.cpp
@@ -103,9 +103,12 @@ std::vector<std::vector<std::size_t>> identifySubdomainMeshElements(
             bulk_mesh.getElementsConnectedToNode(node_id) |
             MeshLib::views::ids | ranges::to<std::vector>;
     }
+
+    auto const& elements = subdomain_mesh.getElements();
 #pragma omp parallel for
-    for (auto* const e : subdomain_mesh.getElements())
+    for (std::ptrdiff_t j = 0; j < std::ssize(elements); ++j)
     {
+        auto* const e = elements[j];
         std::vector<std::size_t> element_node_ids(e->getNumberOfBaseNodes());
         for (unsigned n = 0; n < e->getNumberOfBaseNodes(); ++n)
         {
diff --git a/MeshGeoToolsLib/MeshNodeSearcher.cpp b/MeshGeoToolsLib/MeshNodeSearcher.cpp
index 75b3f1e69ca..a8fdda23c5d 100644
--- a/MeshGeoToolsLib/MeshNodeSearcher.cpp
+++ b/MeshGeoToolsLib/MeshNodeSearcher.cpp
@@ -135,12 +135,13 @@ std::vector<std::size_t> MeshNodeSearcher::getMeshNodeIDs(
 {
     double const epsilon_radius = _search_length_algorithm->getSearchLength();
 
-    std::vector<std::size_t> node_ids;
-    node_ids.reserve(points.size());
+    std::vector<std::size_t> node_ids(points.size());
 
-    for (auto const* const p_ptr : points)
+    auto const number_of_points = std::ssize(points);
+#pragma omp for
+    for (std::ptrdiff_t i = 0; i < number_of_points; ++i)
     {
-        auto const& p = *p_ptr;
+        auto const& p = *points[i];
         std::vector<std::size_t> const ids =
             _mesh_grid.getPointsInEpsilonEnvironment(p, epsilon_radius);
         if (ids.empty())
@@ -168,7 +169,7 @@ std::vector<std::size_t> MeshNodeSearcher::getMeshNodeIDs(
                 ids.size(), p.getID(), p[0], p[1], p[2], epsilon_radius,
                 _mesh.getName(), ss.str());
         }
-        node_ids.push_back(ids.front());
+        node_ids[i] = ids.front();
     }
     return node_ids;
 }
-- 
GitLab