Merge pull request #1904 from endJunction/FindElementsInRadius

Find elements in radius algorithm

MeshLib/findElementsWithinRadius.cpp

+/**
+ * \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)};
+}
+}

MeshLib/findElementsWithinRadius.h

+/**
+ * \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

Tests/MeshLib/TestFindElementsWithinRadius.cpp

+/**
+ * @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);
+}