Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : #include "libmesh/libmesh_config.h"
20 :
21 : // This class is not defined unless libmesh is built with Nanoflann support
22 : #ifdef LIBMESH_HAVE_NANOFLANN
23 :
24 : // libmesh includes
25 : #include "libmesh/point_locator_nanoflann.h"
26 : #include "libmesh/elem.h"
27 : #include "libmesh/libmesh_logging.h"
28 : #include "libmesh/mesh_base.h"
29 : #include "libmesh/mesh_tools.h"
30 : #include "libmesh/int_range.h" // make_range
31 :
32 : // C++ includes
33 : #include <array>
34 : #include <numeric> // std::iota
35 : #include <algorithm> // std::sort
36 :
37 : namespace libMesh
38 : {
39 :
40 0 : PointLocatorNanoflann::PointLocatorNanoflann (const MeshBase & mesh,
41 0 : const PointLocatorBase * master) :
42 : PointLocatorBase (mesh, master),
43 0 : _out_of_mesh_mode(false),
44 0 : _num_results(8)
45 : {
46 0 : this->init();
47 0 : }
48 :
49 :
50 :
51 0 : PointLocatorNanoflann::~PointLocatorNanoflann () = default;
52 :
53 :
54 :
55 : void
56 0 : PointLocatorNanoflann::clear ()
57 : {
58 0 : this->_initialized = false;
59 0 : this->_out_of_mesh_mode = false;
60 :
61 : // reset() actually frees the memory if we are master, otherwise it
62 : // just reduces the ref. count.
63 0 : _elems.reset();
64 0 : _point_cloud.reset();
65 0 : _kd_tree.reset();
66 0 : }
67 :
68 :
69 :
70 : void
71 0 : PointLocatorNanoflann::init ()
72 : {
73 0 : LOG_SCOPE("init()", "PointLocatorNanoflann");
74 :
75 0 : if (!_initialized)
76 : {
77 : // If _master == nullptr, then we _are_ the master, and thus
78 : // responsible for initializing.
79 0 : bool we_are_master = (_master == nullptr);
80 :
81 : // If we are the master PointLocator, fill in the _point_cloud
82 : // data structure with active, local element vertex averages.
83 0 : if (we_are_master)
84 : {
85 0 : _elems = std::make_shared<std::vector<const Elem *>>();
86 0 : _point_cloud = std::make_shared<std::vector<Point>>();
87 :
88 : // Make the KD-Tree out of active element vertex averages.
89 : //
90 : // Note: we use active elements rather than active+local
91 : // elements since we sometimes need to be able to locate
92 : // points in ghosted elements (for example in Periodic BCs).
93 : // Active elements are also natural to use in the
94 : // DistributedMesh case. The size of the KD-Tree on each
95 : // processor therefore scales with the number of elements on
96 : // each processor in either case.
97 : //
98 : // Note 2: the approximate amount of space we should reserve
99 : // here is going to be different for ReplicatedMesh
100 : // (n_active_elem()) vs. DistributedMesh (n_active_local_elem())
101 : // so we just let the containers resize themselves
102 : // automatically.
103 0 : for (const auto & elem : _mesh.active_element_ptr_range())
104 : {
105 0 : _elems->push_back(elem);
106 0 : _point_cloud->push_back(elem->vertex_average());
107 0 : }
108 :
109 : // Construct the KD-Tree
110 : _kd_tree = std::make_shared<kd_tree_t>
111 0 : (LIBMESH_DIM, *this, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));
112 :
113 0 : _kd_tree->buildIndex();
114 : }
115 : else // we are not master
116 : {
117 : // Cast the _master to the appropriate type, and point to
118 : // its data arrays.
119 : const auto my_master =
120 0 : cast_ptr<const PointLocatorNanoflann *>(this->_master);
121 :
122 : // Point our data structures at the master's
123 0 : _elems = my_master->_elems;
124 0 : _point_cloud = my_master->_point_cloud;
125 0 : _kd_tree = my_master->_kd_tree;
126 : }
127 :
128 : // We are initialized now
129 0 : this->_initialized = true;
130 : }
131 0 : }
132 :
133 : nanoflann::KNNResultSet<Real>
134 0 : PointLocatorNanoflann::kd_tree_find_neighbors(const Point & p,
135 : std::size_t num_results) const
136 : {
137 : // We are searching for the Point(s) closest to Point p.
138 : //
139 : // TODO: The kd_tree's findNeighbors() routine needs a pointer to
140 : // Real of length LIBMESH_DIM. It might be convenient if libMesh
141 : // Points had a data() member that provided this, for now we just
142 : // copy the coordinates into a std::array of the appropriate size.
143 : std::array<Real, LIBMESH_DIM> query_pt;
144 0 : for (int i=0; i<LIBMESH_DIM; ++i)
145 0 : query_pt[i] = p(i);
146 :
147 : // Allocate storage for the indices and distances
148 0 : _ret_index.resize(num_results);
149 0 : _out_dist_sqr.resize(num_results);
150 :
151 : // nanoflann::KNNResultSet cannot be resized/reused easily, I think
152 : // it is just meant to be re-created for each search.
153 0 : nanoflann::KNNResultSet<Real> result_set(num_results);
154 :
155 : // Initialize the result_set
156 0 : result_set.init(_ret_index.data(), _out_dist_sqr.data());
157 :
158 : // Do the search
159 : // We leave all the SearchParams() ctor args on their default values
160 : // (note that the first arg is ignored anyway).
161 0 : _kd_tree->findNeighbors(result_set, query_pt.data(), nanoflann::SearchParameters());
162 :
163 0 : return result_set;
164 : }
165 :
166 :
167 :
168 : const Elem *
169 0 : PointLocatorNanoflann::operator() (const Point & p,
170 : const std::set<subdomain_id_type> * allowed_subdomains) const
171 : {
172 0 : libmesh_assert (this->_initialized);
173 :
174 0 : LOG_SCOPE("operator()", "PointLocatorNanoflann");
175 :
176 : // Keep track of the number of elements checked in detail. This can
177 : // be useful for determining an optimal _num_results for
178 : // specific applications, by looking at the average and best/worse
179 : // case of each linear search.
180 0 : unsigned int n_elems_checked = 0;
181 :
182 : // If a containing Elem is found locally, we will set this pointer.
183 0 : const Elem * found_elem = nullptr;
184 :
185 0 : auto result_set = this->kd_tree_find_neighbors(p, _num_results);
186 :
187 : // The results from Nanoflann are already sorted by (squared)
188 : // distance, but within that list of results, there may be some
189 : // vertex averages which are equidistant from the searched-for
190 : // Point. Therefore, we will now indirect_sort the results based on
191 : // elem id, so that the lowest-id Elem from the class of Elems which
192 : // are the same distance away from the search Point is always
193 : // selected.
194 :
195 : // operator< comparison lambda used in the indirect sort. The
196 : // inputs to this vector are indices of the result_set array.
197 0 : auto comp = [this](std::size_t lhs, std::size_t rhs) -> bool
198 : {
199 : // First we sort by squared distance
200 0 : if (_out_dist_sqr[lhs] < _out_dist_sqr[rhs])
201 0 : return true;
202 0 : if (_out_dist_sqr[rhs] < _out_dist_sqr[lhs])
203 0 : return false;
204 :
205 : // If we made it here without returning, then the Points were
206 : // equidistant, so we sort based on Elem id instead.
207 0 : auto nanoflann_index_lhs = _ret_index[lhs];
208 0 : auto elem_id_lhs = (*_elems)[nanoflann_index_lhs]->id();
209 :
210 0 : auto nanoflann_index_rhs = _ret_index[rhs];
211 0 : auto elem_id_rhs = (*_elems)[nanoflann_index_rhs]->id();
212 :
213 0 : if (elem_id_lhs < elem_id_rhs)
214 0 : return true;
215 :
216 0 : return false;
217 0 : };
218 :
219 : // Set up the indices and do the indirect sort. The results are
220 : // stored in _b so that _b[0] is the first result to check, _b[1]
221 : // is second, etc.
222 0 : _b.resize(result_set.size());
223 0 : std::iota(_b.begin(), _b.end(), 0);
224 0 : std::sort(_b.begin(), _b.end(), comp);
225 :
226 : // Note: we use result_set.size() here rather than
227 : // _num_results, in case the result returns fewer than
228 : // _num_results results!
229 0 : for (auto r : make_range(result_set.size()))
230 : {
231 : // Translate the Nanoflann index, which is from [0..n_points),
232 : // into the corresponding Elem id from the mesh. Note that we
233 : // use index _b[r] instead of r, since we want to process the
234 : // Elems in order based on Elem id.
235 0 : auto nanoflann_index = _ret_index[_b[r]];
236 0 : const Elem * candidate_elem = (*_elems)[nanoflann_index];
237 :
238 : // Debugging: print the results
239 : // libMesh::err << "Vertex average/Elem id = " << candidate_elem->id()
240 : // << ", dist^2 = " << _out_dist_sqr[r]
241 : // << std::endl;
242 :
243 : // Before we even check whether the candidate Elem actually
244 : // contains the Point, we may need to check whether the
245 : // candidate Elem is from an allowed subdomain. If the
246 : // candidate Elem is not from an allowed subdomain, we continue
247 : // to the next one.
248 0 : if (allowed_subdomains && !allowed_subdomains->count(candidate_elem->subdomain_id()))
249 : {
250 : // Debugging
251 : // libMesh::err << "Elem " << candidate_elem->id() << " was not from an allowed subdomain, continuing search." << std::endl;
252 0 : continue;
253 : }
254 :
255 : // If we made it here, then the candidate Elem is from an
256 : // allowed subdomain, so let's next check whether it contains
257 : // the point. If the user set a custom tolerance, then we
258 : // actually check close_to_point() rather than contains_point(),
259 : // since this latter function warns about using non-default
260 : // tolerances, but otherwise does the same test.
261 0 : bool inside = _use_contains_point_tol ?
262 0 : candidate_elem->close_to_point(p, _contains_point_tol) :
263 0 : candidate_elem->contains_point(p);
264 :
265 : // Increment the number of elements checked
266 0 : n_elems_checked++;
267 :
268 : // If the point is inside an Elem from an allowed subdomain, we are done.
269 0 : if (inside)
270 : {
271 : // Debugging: report the number of Elems checked
272 : // libMesh::err << "Checked " << n_elems_checked << " nearby Elems before finding a containing Elem." << std::endl;
273 :
274 0 : found_elem = candidate_elem;
275 0 : break;
276 : }
277 : } // end for(r)
278 :
279 : // If we made it here, then at least one of the following happened:
280 : // .) All the candidate elements were from non-allowed subdomains.
281 : // .) The Point was not inside _any_ of the _num_results candidate Elems.
282 : // Thus, if we are not in _out_of_mesh_mode, throw an error,
283 : // otherwise return nullptr to indicate that no suitable element was
284 : // found.
285 0 : if (!_out_of_mesh_mode && !found_elem)
286 : {
287 : // Debugging: we are about to throw an error, but before we do,
288 : // print information about the closest elements (by vertex average
289 : // distance) that the Point was not found in.
290 : // for (auto r : make_range(result_set.size()))
291 : // {
292 : // auto nanoflann_index = _ret_index[_b[r]];
293 : // const Elem * candidate_elem = (*_elems)[nanoflann_index];
294 : //
295 : // libMesh::err << "Vertex average/Elem id = " << candidate_elem->id()
296 : // << ", dist = " << std::sqrt(_out_dist_sqr[_b[r]])
297 : // << std::endl;
298 : // } // end for(r)
299 :
300 :
301 0 : libmesh_error_msg("Point " << p << " was not contained within the closest " << n_elems_checked <<
302 : " elems (by vertex average distance), and _out_of_mesh_mode was not enabled.");
303 : }
304 :
305 0 : return found_elem;
306 : }
307 :
308 :
309 : void
310 0 : PointLocatorNanoflann::operator() (const Point & p,
311 : std::set<const Elem *> & candidate_elements,
312 : const std::set<subdomain_id_type> * allowed_subdomains) const
313 : {
314 0 : libmesh_assert (this->_initialized);
315 :
316 0 : LOG_SCOPE("operator() returning set", "PointLocatorNanoflann");
317 :
318 : // Make sure the output set is initially empty
319 0 : candidate_elements.clear();
320 :
321 : // Keep track of the number of elements checked in detail
322 : // unsigned int n_elems_checked = 0;
323 :
324 : // Do the KD-Tree search
325 0 : auto result_set = this->kd_tree_find_neighbors(p, _num_results);
326 :
327 0 : for (auto r : make_range(result_set.size()))
328 : {
329 : // Translate the Nanoflann index, which is from [0..n_points),
330 : // into the corresponding Elem id from the mesh.
331 0 : auto nanoflann_index = _ret_index[r];
332 0 : const Elem * candidate_elem = (*_elems)[nanoflann_index];
333 :
334 : // Before we even check whether the candidate Elem actually
335 : // contains the Point, we may need to check whether the
336 : // candidate Elem is from an allowed subdomain. If the
337 : // candidate Elem is not from an allowed subdomain, we continue
338 : // to the next one.
339 0 : if (allowed_subdomains && !allowed_subdomains->count(candidate_elem->subdomain_id()))
340 0 : continue;
341 :
342 : // If we made it here, then the candidate Elem is from an
343 : // allowed subdomain, so let's next check whether it contains
344 : // the point. If the user set a custom tolerance, then we
345 : // actually check close_to_point() rather than contains_point(),
346 : // since this latter function warns about using non-default
347 : // tolerances, but otherwise does the same test.
348 0 : bool inside = _use_contains_point_tol ?
349 0 : candidate_elem->close_to_point(p, _contains_point_tol) :
350 0 : candidate_elem->contains_point(p);
351 :
352 : // Increment the number of elements checked
353 : // n_elems_checked++;
354 :
355 : // If the point is contained in/close to an Elem from an
356 : // allowed subdomain, add it to the list.
357 0 : if (inside)
358 0 : candidate_elements.insert(candidate_elem);
359 : } // end for(r)
360 :
361 : // Debugging: for performance reasons, it may be useful to print the
362 : // number of Elems actually checked during the search.
363 : // libMesh::err << "Checked " << n_elems_checked << " nearby Elems before finding a containing Elem." << std::endl;
364 0 : }
365 :
366 :
367 : void
368 0 : PointLocatorNanoflann::enable_out_of_mesh_mode ()
369 : {
370 : // Out-of-mesh mode should now work properly even on meshes with
371 : // non-affine elements.
372 0 : _out_of_mesh_mode = true;
373 0 : }
374 :
375 :
376 : void
377 0 : PointLocatorNanoflann::disable_out_of_mesh_mode ()
378 : {
379 0 : _out_of_mesh_mode = false;
380 0 : }
381 :
382 :
383 :
384 : std::size_t
385 0 : PointLocatorNanoflann::get_num_results() const
386 : {
387 0 : return _num_results;
388 : }
389 :
390 :
391 :
392 : void
393 0 : PointLocatorNanoflann::set_num_results(std::size_t val)
394 : {
395 : // Must request at least 1 result
396 0 : _num_results = std::max(static_cast<std::size_t>(1), val);
397 0 : }
398 :
399 : //
400 : // Required Nanoflann APIs
401 : //
402 :
403 0 : std::size_t PointLocatorNanoflann::kdtree_get_point_count() const
404 : {
405 0 : return _point_cloud->size();
406 : }
407 :
408 :
409 :
410 : PointLocatorNanoflann::coord_t
411 0 : PointLocatorNanoflann::kdtree_distance(const coord_t * p1,
412 : const std::size_t idx_p2,
413 : std::size_t libmesh_dbg_var(size)) const
414 : {
415 : // We only consider LIBMESH_DIM dimensional KD-Trees, so just make
416 : // sure we were called consistently.
417 0 : libmesh_assert_equal_to(size, LIBMESH_DIM);
418 :
419 : // Construct a libmesh Point object from the LIBMESH_DIM-dimensional
420 : // input object, p1.
421 0 : Point point1;
422 0 : for (int i=0; i<LIBMESH_DIM; ++i)
423 0 : point1(i) = p1[i];
424 :
425 : // Compute Euclidean distance, squared
426 0 : return (point1 - (*_point_cloud)[idx_p2]).norm_sq();
427 : }
428 :
429 :
430 :
431 : PointLocatorNanoflann::coord_t
432 0 : PointLocatorNanoflann::kdtree_get_pt(const std::size_t idx, int dim) const
433 : {
434 0 : libmesh_assert_less (idx, _point_cloud->size());
435 0 : libmesh_assert_less (dim, LIBMESH_DIM);
436 :
437 0 : return (*_point_cloud)[idx](dim);
438 : }
439 :
440 : } // namespace libMesh
441 :
442 : #endif
|
