diff --git a/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp
index 71f2575145ce3e14ddf8e2b3ece0d8c9a0cd4af0..91a44c2031a4ca1babae1fc91b4f0e43da6525e5 100644
--- a/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp
+++ b/MeshGeoToolsLib/BoundaryElementsAtPoint.cpp
@@ -9,6 +9,10 @@
 
 #include "BoundaryElementsAtPoint.h"
 
+#ifdef USE_PETSC
+#include <mpi.h>
+#endif
+
 #include "GeoLib/Point.h"
 #include "MathLib/Point3d.h"
 #include "MeshGeoToolsLib/MeshNodeSearcher.h"
@@ -25,6 +29,27 @@ BoundaryElementsAtPoint::BoundaryElementsAtPoint(
     : _mesh(mesh), _point(point)
 {
     auto const node_ids = mshNodeSearcher.getMeshNodeIDs(_point);
+
+#ifdef USE_PETSC
+    std::size_t const number_of_found_nodes_at_rank = node_ids.size();
+    std::size_t number_of_total_found_nodes = 0;
+
+    MPI_Allreduce(&number_of_found_nodes_at_rank, &number_of_total_found_nodes,
+                  1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD);
+
+    if (number_of_total_found_nodes == 0)
+    {
+        OGS_FATAL(
+            "BoundaryElementsAtPoint: the mesh node searcher was unable to "
+            "locate the point ({:f}, {:f}, {:f}) in the mesh.",
+            _point[0], _point[1], _point[2]);
+    }
+
+    if (number_of_found_nodes_at_rank == 0)
+    {
+        return;
+    }
+#else
     if (node_ids.empty())
     {
         OGS_FATAL(
@@ -32,10 +57,13 @@ BoundaryElementsAtPoint::BoundaryElementsAtPoint(
             "locate the point ({:f}, {:f}, {:f}) in the mesh.",
             _point[0], _point[1], _point[2]);
     }
+#endif
+
     if (node_ids.size() == 1)
     {
         std::array<MeshLib::Node*, 1> const nodes = {
             {const_cast<MeshLib::Node*>(_mesh.getNode(node_ids[0]))}};
+
         _boundary_elements.push_back(new MeshLib::Point{nodes, node_ids[0]});
         return;
     }
diff --git a/MeshGeoToolsLib/ConstructMeshesFromGeometries.cpp b/MeshGeoToolsLib/ConstructMeshesFromGeometries.cpp
index 2cd477582bd47db1159e1918b4b5f37f7ca488ee..96b509e5c82d4bb9e8fc71577585ca072bb189e8 100644
--- a/MeshGeoToolsLib/ConstructMeshesFromGeometries.cpp
+++ b/MeshGeoToolsLib/ConstructMeshesFromGeometries.cpp
@@ -9,6 +9,10 @@
 
 #include "ConstructMeshesFromGeometries.h"
 
+#ifdef USE_PETSC
+#include <mpi.h>
+#endif
+
 #include "BaseLib/Logging.h"
 #include "BoundaryElementsSearcher.h"
 #include "GeoLib/GEOObjects.h"
@@ -65,6 +69,10 @@ constructAdditionalMeshesFromGeometries(
                 MeshLib::cloneElements(
                     boundary_element_searcher.getBoundaryElements(
                         geometry, multiple_nodes_allowed))));
+
+#ifdef USE_PETSC
+            MPI_Barrier(MPI_COMM_WORLD);
+#endif
         }
     }
     return additional_meshes;
diff --git a/Tests/MeshGeoToolsLib/TestConstructAdditionalMeshesFromGeoObjects.cpp b/Tests/MeshGeoToolsLib/TestConstructAdditionalMeshesFromGeoObjects.cpp
index 7b387d7277dd56bf0d3d151fbe52181da61d9d8c..14c5b5a4764afac47789679fad845fb0115bcfd9 100644
--- a/Tests/MeshGeoToolsLib/TestConstructAdditionalMeshesFromGeoObjects.cpp
+++ b/Tests/MeshGeoToolsLib/TestConstructAdditionalMeshesFromGeoObjects.cpp
@@ -17,6 +17,11 @@
 #include "MeshGeoToolsLib/ConstructMeshesFromGeometries.h"
 #include "MeshGeoToolsLib/SearchLength.h"
 #include "MeshLib/Mesh.h"
+
+#ifdef USE_PETSC
+#include "MeshLib/NodePartitionedMesh.h"
+#endif
+
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
 #include "Tests/GeoLib/CreateTestPoints.h"
 
@@ -25,8 +30,16 @@ TEST(ConstructAdditionalMeshesFromGeoObjects, PointMesh)
     // create 10x10x10 mesh using 1000 hexahedra
     const double length = 10.0;
     const std::size_t n_subdivisions = 10;
+
+#ifdef USE_PETSC
+    std::unique_ptr<MeshLib::NodePartitionedMesh> mesh =
+        std::make_unique<MeshLib::NodePartitionedMesh>(
+            *MeshLib::MeshGenerator::generateRegularHexMesh(length,
+                                                            n_subdivisions));
+#else
     std::unique_ptr<MeshLib::Mesh> mesh(
         MeshLib::MeshGenerator::generateRegularHexMesh(length, n_subdivisions));
+#endif
 
     // create geometry: for every mesh node exactly one point
     GeoLib::GEOObjects geometries;
@@ -55,8 +68,16 @@ TEST(ConstructAdditionalMeshesFromGeoObjects, PointMeshLargeSearchRadius)
     // create 10x10x10 mesh using 1000 hexahedra
     const double length = 10.0;
     const std::size_t n_subdivisions = 10;
+#ifdef USE_PETSC
+
+    std::unique_ptr<MeshLib::NodePartitionedMesh> mesh =
+        std::make_unique<MeshLib::NodePartitionedMesh>(
+            *MeshLib::MeshGenerator::generateRegularHexMesh(length,
+                                                            n_subdivisions));
+#else
     std::unique_ptr<MeshLib::Mesh> mesh(
         MeshLib::MeshGenerator::generateRegularHexMesh(length, n_subdivisions));
+#endif
 
     // create geometry: for every mesh node exactly one point
     GeoLib::GEOObjects geometries;