preCICE v3.1.2
Loading...
Searching...
No Matches
Index.cpp
Go to the documentation of this file.
1#include <Eigen/Core>
2#include <algorithm>
3#include <boost/iterator/function_output_iterator.hpp>
4#include <boost/range/irange.hpp>
5#include <utility>
6
9#include "profiling/Event.hpp"
10#include "query/Index.hpp"
12
13namespace precice::query {
14
16
17namespace bg = boost::geometry;
18namespace bgi = boost::geometry::index;
19
24
31
44
46{
47 if (indices.vertexRTree) {
48 return indices.vertexRTree;
49 }
50
51 precice::profiling::Event e("query.index.getVertexIndexTree." + mesh.getName());
52
53 // Generating the rtree is expensive, so passing everything in the ctor is
54 // the best we can do. Even passing an index range instead of calling
55 // tree->insert repeatedly is about 10x faster.
58 auto tree = std::make_shared<VertexTraits::RTree>(
59 boost::irange<std::size_t>(0lu, mesh.nVertices()), params, ind);
60
61 indices.vertexRTree = std::move(tree);
62 return indices.vertexRTree;
63}
64
66{
67 if (indices.edgeRTree) {
68 return indices.edgeRTree;
69 }
70
71 precice::profiling::Event e("query.index.getEdgeIndexTree." + mesh.getName());
72
73 // Generating the rtree is expensive, so passing everything in the ctor is
74 // the best we can do. Even passing an index range instead of calling
75 // tree->insert repeatedly is about 10x faster.
78
79 auto tree = std::make_shared<EdgeTraits::RTree>(
80 boost::irange<std::size_t>(0lu, mesh.edges().size()), params, ind);
81
82 indices.edgeRTree = std::move(tree);
83 return indices.edgeRTree;
84}
85
87{
88 if (indices.triangleRTree) {
89 return indices.triangleRTree;
90 }
91
92 precice::profiling::Event e("query.index.getTriangleIndexTree." + mesh.getName());
93
94 // We first generate the values for the triangle rtree.
95 // The resulting vector is a random access range, which can be passed to the
96 // constructor of the rtree for more efficient indexing.
98 elements.reserve(mesh.triangles().size());
99 for (size_t i = 0; i < mesh.triangles().size(); ++i) {
100 auto box = bg::return_envelope<RTreeBox>(mesh.triangles()[i]);
101 elements.emplace_back(std::move(box), i);
102 }
103
104 // Generating the rtree is expensive, so passing everything in the ctor is
105 // the best we can do.
108
109 auto tree = std::make_shared<TriangleTraits::RTree>(elements, params, ind);
110 indices.triangleRTree = std::move(tree);
111 return indices.triangleRTree;
112}
113
115{
116 if (indices.tetraRTree) {
117 return indices.tetraRTree;
118 }
119
120 precice::profiling::Event e("query.index.getTetraIndexTree." + mesh.getName());
121
122 // We first generate the values for the tetra rtree.
123 // The resulting vector is a random access range, which can be passed to the
124 // constructor of the rtree for more efficient indexing.
126 elements.reserve(mesh.tetrahedra().size());
127 for (size_t i = 0; i < mesh.tetrahedra().size(); ++i) {
128 // We use a custom function to compute the AABB, because
129 // bg::return_envelope was designed for polygons.
130 auto box = makeBox(mesh.tetrahedra()[i]);
131 elements.emplace_back(std::move(box), i);
132 }
133
134 // Generating the rtree is expensive, so passing everything in the ctor is
135 // the best we can do.
138
139 auto tree = std::make_shared<TetrahedronTraits::RTree>(elements, params, ind);
140 indices.tetraRTree = std::move(tree);
141 return indices.tetraRTree;
142}
143
145{
146 indices.vertexRTree.reset();
147 indices.edgeRTree.reset();
148 indices.triangleRTree.reset();
149 indices.tetraRTree.reset();
150}
151
152//
153// query::Index
154//
155
157 : _mesh(mesh.get())
158{
159 _pimpl = std::make_unique<IndexImpl>(IndexImpl{});
160}
161
163 : _mesh(&mesh)
164{
165 _pimpl = std::make_unique<IndexImpl>(IndexImpl{});
166}
167
168// Required for the pimpl idiom to work with std::unique_ptr
169Index::~Index() = default;
170
171VertexMatch Index::getClosestVertex(const Eigen::VectorXd &sourceCoord)
172{
174
176 VertexMatch match;
177 const auto &rtree = _pimpl->getVertexRTree(*_mesh);
178 rtree->query(bgi::nearest(sourceCoord, 1), boost::make_function_output_iterator([&](size_t matchID) {
179 match = VertexMatch(matchID);
180 }));
181 return match;
182}
183
184std::vector<VertexID> Index::getClosestVertices(const Eigen::VectorXd &sourceCoord, int n)
185{
188 std::vector<VertexID> matches;
189 const auto & rtree = _pimpl->getVertexRTree(*_mesh);
190
191 rtree->query(bgi::nearest(sourceCoord, n), boost::make_function_output_iterator([&](size_t matchID) {
192 matches.emplace_back(matchID);
193 }));
194 return matches;
195}
196
197std::vector<EdgeMatch> Index::getClosestEdges(const Eigen::VectorXd &sourceCoord, int n)
198{
200
201 const auto &rtree = _pimpl->getEdgeRTree(*_mesh);
202
204 rtree->query(bgi::nearest(sourceCoord, n), boost::make_function_output_iterator([&](size_t matchID) {
205 matches.emplace_back(matchID);
206 }));
207 return matches;
208}
209
210std::vector<TriangleMatch> Index::getClosestTriangles(const Eigen::VectorXd &sourceCoord, int n)
211{
213 const auto &rtree = _pimpl->getTriangleRTree(*_mesh);
214
216 rtree->query(bgi::nearest(sourceCoord, n),
217 boost::make_function_output_iterator([&](TriangleTraits::IndexType const &match) {
218 matches.emplace_back(match.second);
219 }));
220 return matches;
221}
222
224{
226
227 // Prepare boost::geometry box
228 auto coords = centerVertex.getCoords();
229 auto searchBox = query::makeBox(coords.array() - radius, coords.array() + radius);
230
231 const auto & rtree = _pimpl->getVertexRTree(*_mesh);
232 std::vector<VertexID> matches;
233 rtree->query(bgi::intersects(searchBox) and bg::index::satisfies([&](size_t const i) { return bg::distance(centerVertex, _mesh->vertex(i)) < radius; }),
234 std::back_inserter(matches));
235 return matches;
236}
237
238bool Index::isAnyVertexInsideBox(const mesh::Vertex &centerVertex, double radius)
239{
241
242 // Prepare boost::geometry box
243 auto coords = centerVertex.getCoords();
244 auto searchBox = query::makeBox(coords.array() - radius, coords.array() + radius);
245
246 const auto &rtree = _pimpl->getVertexRTree(*_mesh);
247
248 auto queryIter = rtree->qbegin(bgi::intersects(searchBox) and bg::index::satisfies([&](size_t const i) { return bg::distance(centerVertex, _mesh->vertex(i)) < radius; }));
249 bool hasMatch = queryIter != rtree->qend();
250 return hasMatch;
251}
252
254{
256 // Add tree to the local cache
257 const auto & rtree = _pimpl->getVertexRTree(*_mesh);
258 std::vector<VertexID> matches;
259 rtree->query(bgi::intersects(query::makeBox(bb.minCorner(), bb.maxCorner())), std::back_inserter(matches));
260 return matches;
261}
262
264{
266 const auto &rtree = _pimpl->getTetraRTree(*_mesh);
267
269 rtree->query(bgi::covers(location), boost::make_function_output_iterator([&](TetrahedronTraits::IndexType const &match) {
270 matches.emplace_back(match.second);
271 }));
272 return matches;
273}
274
275ProjectionMatch Index::findNearestProjection(const Eigen::VectorXd &location, int n)
276{
277 if (_mesh->getDimensions() == 2) {
278 return findEdgeProjection(location, n, findVertexProjection(location));
279 } else {
280 return findTriangleProjection(location, n, findVertexProjection(location));
281 }
282}
283
284ProjectionMatch Index::findCellOrProjection(const Eigen::VectorXd &location, int n)
285{
286 if (_mesh->getDimensions() == 2) {
287 auto matchedTriangles = getClosestTriangles(location, n);
288 for (const auto &match : matchedTriangles) {
289 auto polation = mapping::Polation(location, _mesh->triangles()[match.index]);
290 if (polation.isInterpolation()) {
291 return {std::move(polation)};
292 }
293 }
294
295 // If no triangle is found, fall-back on NP
296 return findNearestProjection(location, n);
297 } else {
298
299 // Find correct tetra, or fall back to NP
300 auto matchedTetra = getEnclosingTetrahedra(location);
301 for (const auto &match : matchedTetra) {
302 // Matches are raw indices, not (indices, distance) pairs
303 auto polation = mapping::Polation(location, _mesh->tetrahedra()[match]);
304 if (polation.isInterpolation()) {
305 return {std::move(polation)};
306 }
307 }
308 return findNearestProjection(location, n);
309 }
310}
311
312ProjectionMatch Index::findVertexProjection(const Eigen::VectorXd &location)
313{
314 auto match = getClosestVertex(location);
315 return {mapping::Polation{location, _mesh->vertex(match.index)}};
316}
317
318ProjectionMatch Index::findEdgeProjection(const Eigen::VectorXd &location, int n, ProjectionMatch closestVertex)
319{
321 candidates.reserve(n);
322 for (const auto &match : getClosestEdges(location, n)) {
323 auto polation = mapping::Polation(location, _mesh->edges()[match.index]);
324 if (polation.isInterpolation()) {
325 candidates.emplace_back(std::move(polation));
326 }
327 }
328
329 // Could not find edge projection element, fall back to vertex projection
330 if (candidates.empty()) {
331 return closestVertex;
332 }
333
334 // Prefer the closest vertex if it closer than the closest edge
335 auto min = std::min_element(candidates.begin(), candidates.end());
336 if (min->polation.distance() > closestVertex.polation.distance()) {
337 return closestVertex;
338 }
339 return *min;
340}
341
342ProjectionMatch Index::findTriangleProjection(const Eigen::VectorXd &location, int n, ProjectionMatch closestVertex)
343{
345 candidates.reserve(n);
346 for (const auto &match : getClosestTriangles(location, n)) {
347 auto polation = mapping::Polation(location, _mesh->triangles()[match.index]);
348 if (polation.isInterpolation()) {
349 candidates.emplace_back(std::move(polation));
350 }
351 }
352
353 // Could not find triangle projection element, fall back to edge projection
354 if (candidates.empty()) {
355 return findEdgeProjection(location, n, std::move(closestVertex));
356 }
357
358 // Fallback to edge projection if a vertex is closer than the best triangle match
359 auto min = std::min_element(candidates.begin(), candidates.end());
360 if (min->polation.distance() > closestVertex.polation.distance()) {
361 return findEdgeProjection(location, n, std::move(closestVertex));
362 }
363
364 return *min;
365}
366
368{
370 // if the mesh is empty, we will most likely hit an assertion in the bounding box class
371 // therefore, we keep the assert here, but might want to return an empty bounding box in case
372 // we want to allow calling this function with empty meshes
374
375 auto rtreeBox = _pimpl->getVertexRTree(*_mesh)->bounds();
376 int dim = _mesh->getDimensions();
377 Eigen::VectorXd min(dim), max(dim);
378
379 min[0] = rtreeBox.min_corner().get<0>();
380 min[1] = rtreeBox.min_corner().get<1>();
381 max[0] = rtreeBox.max_corner().get<0>();
382 max[1] = rtreeBox.max_corner().get<1>();
383
384 if (dim > 2) {
385 min[2] = rtreeBox.min_corner().get<2>();
386 max[2] = rtreeBox.max_corner().get<2>();
387 }
388 return mesh::BoundingBox{min, max};
389}
390
392{
393 _pimpl->clear();
394}
395
396} // namespace precice::query
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T back_inserter(T... args)
T begin(T... args)
This class provides a lightweight logger.
Definition Logger.hpp:16
Calculates the barycentric coordinates of a coordinate on the given vertex/edge/triangle and stores t...
Definition Polation.hpp:24
double distance() const
Returns the projection distance.
Definition Polation.cpp:91
An axis-aligned bounding box around a (partition of a) mesh.
Eigen::VectorXd maxCorner() const
the max corner of the bounding box
Eigen::VectorXd minCorner() const
the min corner of the bounding box
Container and creator for meshes.
Definition Mesh.hpp:39
int getDimensions() const
Definition Mesh.cpp:98
VertexContainer & vertices()
Returns modifieable container holding all vertices.
Definition Mesh.cpp:53
const std::string & getName() const
Returns the name of the mesh, as set in the config file.
Definition Mesh.cpp:218
std::size_t nVertices() const
Returns the number of vertices.
Definition Mesh.cpp:63
TetraContainer & tetrahedra()
Returns modifiable container holding all tetrahedra.
Definition Mesh.cpp:93
Vertex & vertex(VertexID id)
Mutable access to a vertex by VertexID.
Definition Mesh.cpp:41
bool empty() const
Does the mesh contain any vertices?
Definition Mesh.hpp:88
TriangleContainer & triangles()
Returns modifiable container holding all triangles.
Definition Mesh.cpp:78
EdgeContainer & edges()
Returns modifiable container holding all edges.
Definition Mesh.cpp:68
Vertex of a mesh.
Definition Vertex.hpp:16
Eigen::VectorXd getCoords() const
Returns the coordinates of the vertex.
Definition Vertex.hpp:116
VertexTraits::Ptr getVertexRTree(const mesh::Mesh &mesh)
Definition Index.cpp:45
TetrahedronTraits::Ptr getTetraRTree(const mesh::Mesh &mesh)
Definition Index.cpp:114
TriangleTraits::Ptr getTriangleRTree(const mesh::Mesh &mesh)
Definition Index.cpp:86
EdgeTraits::Ptr getEdgeRTree(const mesh::Mesh &mesh)
Definition Index.cpp:65
mesh::BoundingBox getRtreeBounds()
Definition Index.cpp:367
std::vector< TetrahedronID > getEnclosingTetrahedra(const Eigen::VectorXd &location)
Return all the tetrahedra whose axis-aligned bounding box contains a vertex.
Definition Index.cpp:263
mesh::Mesh * _mesh
The indexed Mesh.
Definition Index.hpp:121
ProjectionMatch findCellOrProjection(const Eigen::VectorXd &location, int n)
Definition Index.cpp:284
VertexMatch getClosestVertex(const Eigen::VectorXd &sourceCoord)
Get the closest vertex to the given vertex.
Definition Index.cpp:171
ProjectionMatch findNearestProjection(const Eigen::VectorXd &location, int n)
Find the closest interpolation element to the given location. If exists, triangle or edge projection ...
Definition Index.cpp:275
bool isAnyVertexInsideBox(const mesh::Vertex &centerVertex, double radius)
Returns.
Definition Index.cpp:238
std::unique_ptr< IndexImpl > _pimpl
Definition Index.hpp:118
std::vector< TriangleMatch > getClosestTriangles(const Eigen::VectorXd &sourceCoord, int n)
Get n number of closest triangles to the given vertex.
Definition Index.cpp:210
std::vector< VertexID > getVerticesInsideBox(const mesh::Vertex &centerVertex, double radius)
Return all the vertices inside the box formed by vertex and radius (boundary exclusive)
Definition Index.cpp:223
void clear()
Clear the index.
Definition Index.cpp:391
ProjectionMatch findEdgeProjection(const Eigen::VectorXd &location, int n, ProjectionMatch closestVertex)
Find closest edge interpolation element. If cannot be found, it falls back to vertex projection.
Definition Index.cpp:318
ProjectionMatch findVertexProjection(const Eigen::VectorXd &location)
Closest vertex projection element is always the nearest neighbor.
Definition Index.cpp:312
Index(mesh::PtrMesh mesh)
Definition Index.cpp:156
std::vector< VertexID > getClosestVertices(const Eigen::VectorXd &sourceCoord, int n)
Get n number of closest vertices to the given vertex.
Definition Index.cpp:184
ProjectionMatch findTriangleProjection(const Eigen::VectorXd &location, int n, ProjectionMatch closestVertex)
Find closest face interpolation element. If cannot be found, it falls back to first edge interpolatio...
Definition Index.cpp:342
std::vector< EdgeMatch > getClosestEdges(const Eigen::VectorXd &sourceCoord, int n)
Get n number of closest edges to the given vertex.
Definition Index.cpp:197
static precice::logging::Logger _log
Definition Index.hpp:123
T emplace_back(T... args)
T empty(T... args)
T end(T... args)
T min_element(T... args)
boost::geometry::index::rstar< 16 > RTreeParameters
The general rtree parameter type used in precice.
contains geometrical queries.
RTreeBox makeBox(const pm::Vertex::RawCoords &min, const pm::Vertex::RawCoords &max)
MatchType< struct VertexMatchTag > VertexMatch
Definition Index.hpp:39
constexpr auto get(span< E, S > s) -> decltype(s[N])
Definition span.hpp:602
T reserve(T... args)
T size(T... args)
Struct to hold the index of a primitive match.
Definition Index.hpp:30
VertexTraits::Ptr vertexRTree
Definition Index.cpp:26
EdgeTraits::Ptr edgeRTree
Definition Index.cpp:27
TetrahedronTraits::Ptr tetraRTree
Definition Index.cpp:29
TriangleTraits::Ptr triangleRTree
Definition Index.cpp:28
Struct representing a projection match.
Definition Index.hpp:45
mapping::Polation polation
Definition Index.hpp:56
The type traits of a rtree based on a Primitive.
typename std::conditional< IsDirectIndexable< Primitive >::value, MeshContainerIndex, std::pair< RTreeBox, MeshContainerIndex > >::type IndexType
typename std::conditional< IsDirectIndexable< Primitive >::value, impl::VectorIndexable< MeshContainer >, boost::geometry::index::indexable< IndexType > >::type IndexGetter