From 4538b689cf672d0b6148e63de7b4108de66e4084 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Tue, 8 Aug 2017 00:16:33 +0200
Subject: [PATCH] [MeL] Implement findElementsInRadius algorithm.
Randomized tests on 2D and 3D meshes are provided.
---
MeshLib/ElementAlgorithms.cpp | 88 ++++++++
MeshLib/ElementAlgorithms.h | 36 ++++
Tests/MeshLib/TestFindElementsInRadius.cpp | 226 +++++++++++++++++++++
3 files changed, 350 insertions(+)
create mode 100644 MeshLib/ElementAlgorithms.cpp
create mode 100644 MeshLib/ElementAlgorithms.h
create mode 100644 Tests/MeshLib/TestFindElementsInRadius.cpp
diff --git a/MeshLib/ElementAlgorithms.cpp b/MeshLib/ElementAlgorithms.cpp
new file mode 100644
index 00000000000..739da0a99b0
--- /dev/null
+++ b/MeshLib/ElementAlgorithms.cpp
@@ -0,0 +1,88 @@
+/**
+ * \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 "ElementAlgorithms.h"
+
+
+#include <deque>
+#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> findElementsInRadius(Element const& start_element,
+ double const radius)
+{
+ // Special case for 0 radius. All other radii will include at least one
+ // neighbor of the start element.
+ if (radius == 0.)
+ return {start_element.getID()};
+
+ // Returns true if the test node is inside the radius of any of the
+ // element's nodes.
+ auto node_inside_radius = [&start_element, radius](Node const* test_node) {
+ for (unsigned n = 0; n < start_element.getNumberOfNodes(); ++n)
+ {
+ if (MathLib::sqrDist(*test_node, *start_element.getNode(n)) <=
+ radius * radius)
+ return true;
+ }
+ return false;
+ };
+
+ std::set<std::size_t> found_elements;
+ std::deque<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());
+ for (unsigned n = 0; n < element.getNumberOfNeighbors(); ++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.front();
+ neighbors_to_visit.pop_front();
+
+ // if any node is inside the radius, all neighbors are visited.
+ for (unsigned i = 0; i < current_element.getNumberOfNodes(); ++i)
+ {
+ if (!node_inside_radius(current_element.getNode(i)))
+ continue;
+
+ mark_visited_and_add_neighbors(current_element);
+ }
+ }
+
+ return {std::begin(found_elements), std::end(found_elements)};
+}
+}
diff --git a/MeshLib/ElementAlgorithms.h b/MeshLib/ElementAlgorithms.h
new file mode 100644
index 00000000000..f92be6651c2
--- /dev/null
+++ b/MeshLib/ElementAlgorithms.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> findElementsInRadius(Element const& e,
+ double const radius);
+
+} // namespace MeshLib
diff --git a/Tests/MeshLib/TestFindElementsInRadius.cpp b/Tests/MeshLib/TestFindElementsInRadius.cpp
new file mode 100644
index 00000000000..4cd4207950d
--- /dev/null
+++ b/Tests/MeshLib/TestFindElementsInRadius.cpp
@@ -0,0 +1,226 @@
+/**
+ * @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/ElementAlgorithms.h"
+#include "MeshLib/Elements/Element.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
+#include "MeshLib/Node.h"
+
+#include "Tests/AutoCheckTools.h"
+
+using namespace MeshLib;
+namespace ac = autocheck;
+
+struct MeshLibFindElementInRadius : public ::testing::Test
+{
+ ac::gtest_reporter gtest_reporter;
+};
+
+// For zero radius only the starting element is alway returned.
+TEST_F(MeshLibFindElementInRadius, ZeroRadius)
+{
+ auto mesh = MeshGenerator::generateRegularQuadMesh(10., 10);
+
+ auto same_element_returned = [&mesh](std::size_t& element_id) -> bool {
+ auto result = findElementsInRadius(*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> findByNodeConnectedElements(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*> findNodesInRadius(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*> findNodesInRadiusOfElement(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 = findNodesInRadius(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> bruteForceFindElementIdsInRadius(
+ Mesh const& mesh, Element const& element, double const radius)
+{
+ std::vector<std::size_t> connected_elements;
+ auto const nodes_in_radius =
+ findNodesInRadiusOfElement(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(MeshLibFindElementInRadius, 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 = findElementsInRadius(element, 1e-5);
+
+ std::sort(std::begin(result), std::end(result));
+ auto const expected_elements = findByNodeConnectedElements(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(MeshLibFindElementInRadius, VeryLargeRadius)
+{
+ auto mesh = *MeshGenerator::generateRegularQuadMesh(10., 10);
+
+ auto all_elements_returned = [&mesh](std::size_t& element_id) -> bool {
+ auto result = findElementsInRadius(*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);
+}
+
+// Random test (element_id and radius > 0); comparison with brute-force search
+// algorithm.
+TEST_F(MeshLibFindElementInRadius, RandomPositiveRadius2d)
+{
+ auto mesh = *MeshGenerator::generateRegularQuadMesh(10., 10);
+
+ auto compare_to_brute_force_search = [&mesh](std::size_t& element_id,
+ double const radius) -> bool {
+ auto const& element = *mesh.getElement(element_id);
+
+ auto result = findElementsInRadius(element, radius);
+ std::sort(std::begin(result), std::end(result));
+
+ auto const expected_elements =
+ bruteForceFindElementIdsInRadius(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));
+ };
+
+ 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(MeshLibFindElementInRadius, RandomPositiveRadius3d)
+{
+ auto mesh = *MeshGenerator::generateRegularHexMesh(10., 10);
+
+ auto compare_to_brute_force_search = [&mesh](std::size_t& element_id,
+ double const radius) -> bool {
+ auto const& element = *mesh.getElement(element_id);
+
+ auto result = findElementsInRadius(element, radius);
+ std::sort(std::begin(result), std::end(result));
+
+ auto const expected_elements =
+ bruteForceFindElementIdsInRadius(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));
+ };
+
+ 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);
+}
--
GitLab