diff --git a/MeshLib/findElementsWithinRadius.cpp b/MeshLib/findElementsWithinRadius.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bc69f4875ac23a76707ef50100cffb8baa477ae9
--- /dev/null
+++ b/MeshLib/findElementsWithinRadius.cpp
@@ -0,0 +1,105 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "findElementsWithinRadius.h"
+
+
+#include <vector>
+#include <unordered_set>
+
+#include "BaseLib/uniqueInsert.h"
+#include "MathLib/MathTools.h"
+
+#include "Elements/Element.h"
+#include "Node.h"
+
+namespace MeshLib {
+std::vector<std::size_t> findElementsWithinRadius(Element const& start_element,
+ double const radius_squared)
+{
+ // Special case for 0 radius. All other radii will include at least one
+ // neighbor of the start element.
+ if (radius_squared == 0.)
+ return {start_element.getID()};
+
+ // Collect start element node coordinates into local contigiuos memory.
+ std::vector<MathLib::Point3d> start_element_nodes;
+ {
+ auto const start_element_n_nodes = start_element.getNumberOfNodes();
+ start_element_nodes.reserve(start_element_n_nodes);
+ for (unsigned n = 0; n < start_element_n_nodes; ++n)
+ start_element_nodes.push_back(*start_element.getNode(n));
+ }
+
+ // Returns true if the test node is inside the radius of any of the
+ // element's nodes.
+ auto node_inside_radius = [&start_element_nodes,
+ radius_squared](Node const* test_node) {
+ for (auto const& n : start_element_nodes)
+ {
+ if (MathLib::sqrDist(*test_node, n) <= radius_squared)
+ return true;
+ }
+ return false;
+ };
+
+ // Returns true if any of the test element's nodes is inside the start
+ // element's radius.
+ auto element_in_radius = [&node_inside_radius](Element const& element) {
+ auto const n_nodes = element.getNumberOfNodes();
+ for (unsigned i = 0; i < n_nodes; ++i)
+ {
+ if (node_inside_radius(element.getNode(i)))
+ return true;
+ }
+ return false;
+ };
+
+ std::set<std::size_t> found_elements;
+ std::vector<Element const*> neighbors_to_visit;
+ std::unordered_set<std::size_t> visited_elements;
+
+ auto mark_visited_and_add_neighbors =
+ [&found_elements, &neighbors_to_visit, &visited_elements](
+ Element const& element) {
+ found_elements.insert(element.getID());
+ visited_elements.insert(element.getID());
+ auto const n_neighbors = element.getNumberOfNeighbors();
+ for (unsigned n = 0; n < n_neighbors; ++n)
+ {
+ auto neighbor = element.getNeighbor(n);
+ if (neighbor == nullptr)
+ continue;
+ auto x = visited_elements.find(neighbor->getID());
+ if (x != visited_elements.end())
+ continue;
+
+ neighbors_to_visit.push_back(neighbor);
+
+ visited_elements.insert(neighbor->getID());
+ }
+ };
+
+ // The result always contains the starting element.
+ mark_visited_and_add_neighbors(start_element);
+
+ while (!neighbors_to_visit.empty())
+ {
+ auto const& current_element = *neighbors_to_visit.back();
+ neighbors_to_visit.pop_back();
+
+ // If any node is inside the radius, all neighbors are visited.
+ if (element_in_radius(current_element))
+ mark_visited_and_add_neighbors(current_element);
+ }
+
+ return {std::begin(found_elements), std::end(found_elements)};
+}
+}
diff --git a/MeshLib/findElementsWithinRadius.h b/MeshLib/findElementsWithinRadius.h
new file mode 100644
index 0000000000000000000000000000000000000000..603dc31273c1156ddf24da0fd168d271e7a14d45
--- /dev/null
+++ b/MeshLib/findElementsWithinRadius.h
@@ -0,0 +1,36 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, 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>
+
+namespace MeshLib
+{
+class Element;
+} // namespace MeshLib
+
+namespace MeshLib
+{
+/// Find neighbor elements in given radius of the element.
+/// \returns ids of elements for which at least one node lies inside the radius
+/// of any node of the given element.
+/// \note There are corner-cases not correctly handled by this algorithm:
+/// Elements where an edge, but not the edge' nodes, fall inside the search
+/// radius are not included. The algorithm is only considering node to node
+/// distances and does not take into account if there is a part of an element
+/// (an edge or a face) but not the nodes falls inside the search radius.
+///
+/// \note For radius 0 only the given element's id is returned.
+std::vector<std::size_t> findElementsWithinRadius(Element const& e,
+ double const radius_squared);
+
+} // namespace MeshLib
diff --git a/Tests/MeshLib/TestFindElementsWithinRadius.cpp b/Tests/MeshLib/TestFindElementsWithinRadius.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dc74065d6aea97a40a5e72698708658284a26379
--- /dev/null
+++ b/Tests/MeshLib/TestFindElementsWithinRadius.cpp
@@ -0,0 +1,224 @@
+/**
+ * @copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/LICENSE.txt
+ */
+
+#include <gtest/gtest.h>
+#include <autocheck/autocheck.hpp>
+
+#include "BaseLib/makeVectorUnique.h"
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
+#include "MeshLib/Node.h"
+#include "MeshLib/findElementsWithinRadius.h"
+
+#include "Tests/AutoCheckTools.h"
+
+using namespace MeshLib;
+namespace ac = autocheck;
+
+struct MeshLibFindElementWithinRadius : public ::testing::Test
+{
+ ac::gtest_reporter gtest_reporter;
+};
+
+// For zero radius only the starting element is alway returned.
+TEST_F(MeshLibFindElementWithinRadius, ZeroRadius)
+{
+ auto mesh = MeshGenerator::generateRegularQuadMesh(10., 10);
+
+ auto same_element_returned = [&mesh](std::size_t& element_id) -> bool {
+ auto result =
+ findElementsWithinRadius(*mesh->getElement(element_id), 0.);
+ return (result.size() == 1) && (result[0] == element_id);
+ };
+
+ ac::check<std::size_t>(same_element_returned, 100,
+ ac::make_arbitrary<std::size_t>(), gtest_reporter);
+}
+
+std::vector<std::size_t> findElementIdsConnectedToNode(Node const& node)
+{
+ std::vector<std::size_t> connected_elements;
+ auto const& elements_of_node = node.getElements();
+ std::transform(std::begin(elements_of_node),
+ std::end(elements_of_node),
+ std::back_inserter(connected_elements),
+ [](auto const& e) { return e->getID(); });
+ return connected_elements;
+}
+
+// Find all elements connected through all nodes of the element
+std::vector<std::size_t> findNextNeighborsConnectedByNode(
+ Element const& element)
+{
+ std::vector<std::size_t> connected_elements;
+ for (unsigned n = 0; n < element.getNumberOfBaseNodes(); ++n)
+ {
+ auto const& node = *element.getNode(n);
+ auto const& elements = findElementIdsConnectedToNode(node);
+ std::copy(std::begin(elements), std::end(elements),
+ std::back_inserter(connected_elements));
+ }
+
+ BaseLib::makeVectorUnique(connected_elements);
+ return connected_elements;
+}
+
+// Find nodes in radius of a given node.
+std::vector<Node const*> findNodesWithinRadius(Mesh const& mesh,
+ Node const& node,
+ double const radius)
+{
+ std::vector<Node const*> nodes;
+
+ for (std::size_t i = 0; i < mesh.getNumberOfNodes(); ++i)
+ {
+ auto const* n = mesh.getNode(i);
+ auto const distance_squared = MathLib::sqrDist(node, *n);
+ if (distance_squared > radius * radius)
+ continue;
+
+ nodes.push_back(n);
+ }
+ return nodes;
+}
+
+// Find nodes in radius of all nodes of a given element.
+std::vector<Node const*> findNodesWithinRadius(Mesh const& mesh,
+ Element const& element,
+ double const radius)
+{
+ std::vector<Node const*> nodes_in_radius;
+
+ for (unsigned n = 0; n < element.getNumberOfBaseNodes(); ++n)
+ {
+ auto const& node = *element.getNode(n);
+ auto const& nodes = findNodesWithinRadius(mesh, node, radius);
+ std::copy(std::begin(nodes), std::end(nodes),
+ std::back_inserter(nodes_in_radius));
+ }
+
+ BaseLib::makeVectorUnique(nodes_in_radius);
+ return nodes_in_radius;
+}
+
+// Brute force search for checking the actual algorithm.
+std::vector<std::size_t> bruteForceFindElementIdsWithinRadius(
+ Mesh const& mesh, Element const& element, double const radius)
+{
+ std::vector<std::size_t> connected_elements;
+ auto const nodes_in_radius = findNodesWithinRadius(mesh, element, radius);
+ for (auto n : nodes_in_radius)
+ {
+ auto const& elements = findElementIdsConnectedToNode(*n);
+ std::copy(std::begin(elements), std::end(elements),
+ std::back_inserter(connected_elements));
+ }
+
+ BaseLib::makeVectorUnique(connected_elements);
+ return connected_elements;
+}
+
+// For a small radius only the element and its neighbors (through all nodes) are
+// expected.
+TEST_F(MeshLibFindElementWithinRadius, VerySmallRadius)
+{
+ auto mesh = MeshGenerator::generateRegularQuadMesh(10., 10);
+
+ auto neighboring_elements_returned =
+ [&mesh](std::size_t& element_id) -> bool {
+ auto const& element = *mesh->getElement(element_id);
+ auto result = findElementsWithinRadius(element, 1e-5);
+
+ std::sort(std::begin(result), std::end(result));
+ auto const expected_elements =
+ findNextNeighborsConnectedByNode(element);
+
+ return result.size() == expected_elements.size() &&
+ std::includes(std::begin(result), std::end(result),
+ std::begin(expected_elements),
+ std::end(expected_elements));
+ };
+
+ ac::check<std::size_t>(neighboring_elements_returned, 100,
+ ac::make_arbitrary<std::size_t>(), gtest_reporter);
+}
+
+// For radii large enough to cover all of the mesh all of the elements are
+// expected to be found.
+TEST_F(MeshLibFindElementWithinRadius, VeryLargeRadius)
+{
+ auto mesh = *MeshGenerator::generateRegularQuadMesh(10., 10);
+
+ auto all_elements_returned = [&mesh](std::size_t& element_id) -> bool {
+ auto result =
+ findElementsWithinRadius(*mesh.getElement(element_id), 1e5);
+ BaseLib::makeVectorUnique(result);
+
+ return result.size() == mesh.getNumberOfElements();
+ };
+
+ ac::check<std::size_t>(all_elements_returned, 100,
+ ac::make_arbitrary<std::size_t>(), gtest_reporter);
+}
+
+// Comparison with brute-force search algorithm.
+struct CompareToBruteForceSearch
+{
+ Mesh const& mesh;
+
+ bool operator()(std::size_t& element_id, double const radius) const
+ {
+ auto const& element = *mesh.getElement(element_id);
+
+ auto result = findElementsWithinRadius(element, radius * radius);
+ std::sort(std::begin(result), std::end(result));
+
+ auto const expected_elements =
+ bruteForceFindElementIdsWithinRadius(mesh, element, radius);
+
+ return result.size() == expected_elements.size() &&
+ std::includes(std::begin(result), std::end(result),
+ std::begin(expected_elements),
+ std::end(expected_elements));
+ }
+};
+
+// Random test (element_id and radius > 0); comparison with brute-force search
+// algorithm.
+TEST_F(MeshLibFindElementWithinRadius, RandomPositiveRadius2d)
+{
+ auto mesh = *MeshGenerator::generateRegularQuadMesh(10., 10);
+
+ auto const compare_to_brute_force_search = CompareToBruteForceSearch{mesh};
+
+ ac::check<std::size_t, double>(
+ compare_to_brute_force_search, 100,
+ ac::make_arbitrary(ac::generator<std::size_t>(),
+ ac::map(&ac::absoluteValue, ac::generator<double>()))
+ .discard_if(
+ [](std::size_t, double const r) { return (r < 1e-16); }),
+ gtest_reporter);
+}
+
+// Random test (element_id and radius > 0); comparison with brute-force search
+// algorithm.
+TEST_F(MeshLibFindElementWithinRadius, RandomPositiveRadius3d)
+{
+ auto mesh = *MeshGenerator::generateRegularHexMesh(10., 10);
+
+ auto const compare_to_brute_force_search = CompareToBruteForceSearch{mesh};
+
+ ac::check<std::size_t, double>(
+ compare_to_brute_force_search, 100,
+ ac::make_arbitrary(ac::generator<std::size_t>(),
+ ac::map(&ac::absoluteValue, ac::generator<double>()))
+ .discard_if(
+ [](std::size_t, double const r) { return (r < 1e-16); }),
+ gtest_reporter);
+}