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 :
20 : #ifndef LIBMESH_POINT_LOCATOR_NANOFLANN_H
21 : #define LIBMESH_POINT_LOCATOR_NANOFLANN_H
22 :
23 : #include "libmesh/libmesh_config.h"
24 :
25 : // This class is not defined unless libmesh is built with Nanoflann support
26 : #ifdef LIBMESH_HAVE_NANOFLANN
27 :
28 : // libmesh includes
29 : #include "libmesh/point_locator_base.h"
30 : #include "libmesh/point.h"
31 : #include "libmesh/bounding_box.h"
32 :
33 : // contrib includes
34 : #include "libmesh/nanoflann.hpp"
35 :
36 : // C++ includes
37 : #include <vector>
38 : #include <memory>
39 :
40 : namespace libMesh
41 : {
42 :
43 : // Forward Declarations
44 : class MeshBase;
45 : class Point;
46 : class Elem;
47 :
48 : /**
49 : * This is a PointLocator that uses Nanoflann for its implementation.
50 : * Nanoflann is distributed with libmesh (see: contrib/nanoflann) and
51 : * libmesh must be built with nanoflann enabled for this class to
52 : * work.
53 : *
54 : * This PointLocator is still considered "experimental", so it is not
55 : * enabled in libMesh by default. In particular, there may be false
56 : * negatives (i.e. failures to find a containing Elem for a given
57 : * Point even though there actually is one in the Mesh) on adaptively
58 : * refined meshes, meshes with intersecting mixed-dimension manifolds,
59 : * and meshes with large aspect ratio changes across neighboring
60 : * elements. That said, the Nanoflann PointLocator did successfully
61 : * pass all the MOOSE CI testing that we threw at it (and even
62 : * uncovered some bugs in said tests) so it will likely work well for
63 : * most applications outside of the categories mentioned above.
64 : *
65 : * You must configure libmesh with --enable-nanoflann-pointlocator in
66 : * order to use this class.
67 : *
68 : * \author John W. Peterson
69 : * \date 2021
70 : */
71 0 : class PointLocatorNanoflann : public PointLocatorBase
72 : {
73 : public:
74 : /**
75 : * Constructor. Needs the \p mesh in which the points should be
76 : * located. Optionally takes a pointer to a "master" PointLocator
77 : * object. If non-nullptr, this object simply forwards its calls
78 : * onto the master, so we can have multiple pointers that use the
79 : * same Nanoflann KD-Tree data structure.
80 : */
81 : PointLocatorNanoflann (const MeshBase & mesh,
82 : const PointLocatorBase * master = nullptr);
83 :
84 : /**
85 : * Destructor.
86 : */
87 : virtual ~PointLocatorNanoflann ();
88 :
89 : /**
90 : * Restore to PointLocator to a just-constructed state.
91 : */
92 : virtual void clear() override final;
93 :
94 : /**
95 : * Initializes the locator, so that the \p operator() methods can
96 : * be used. This function allocates dynamic memory with "new".
97 : */
98 : virtual void init() override final;
99 :
100 : /**
101 : * Locates the element in which the point with global coordinates \p
102 : * p is located, optionally restricted to a set of allowed
103 : * subdomains.
104 : */
105 : virtual const Elem * operator() (const Point & p,
106 : const std::set<subdomain_id_type> * allowed_subdomains = nullptr) const override final;
107 :
108 : /**
109 : * Locates a set of elements in proximity to the point with global
110 : * coordinates \p p. Optionally allows the user to restrict the
111 : * subdomains searched. The idea here is that if a Point lies on
112 : * the boundary between two elements, so that it is "in" both
113 : * elements (to within some tolerance) then the candidate_elements
114 : * set will contain both of these elements instead of just picking
115 : * one or the other.
116 : */
117 : virtual void operator() (const Point & p,
118 : std::set<const Elem *> & candidate_elements,
119 : const std::set<subdomain_id_type> * allowed_subdomains = nullptr) const override final;
120 :
121 : /**
122 : * Enables out-of-mesh mode. In this mode, if a searched-for Point
123 : * is not contained in any element of the Mesh, return nullptr
124 : * instead of throwing an error. By default, this mode is off.
125 : */
126 : virtual void enable_out_of_mesh_mode () override final;
127 :
128 : /**
129 : * Disables out-of-mesh mode (default). See above.
130 : */
131 : virtual void disable_out_of_mesh_mode () override final;
132 :
133 : /**
134 : * Set/get the number of results returned by each Nanoflann
135 : * findNeighbors() search. If a containing Elem (to within
136 : * _contains_point_tol) is not found after linear searching the
137 : * first _num_results Elems returned by Nanoflann, then we
138 : * give up and return nullptr (or throw an error if not in
139 : * out-of-mesh-mode).
140 : *
141 : * Remarks:
142 : * 1.) Although we do a linear search through the Nanoflann results,
143 : * the Nanoflann results are always sorted by distance to the
144 : * searched-for Point, so it's likely that the containing Elem will
145 : * be found at the beginning of the linear search rather than the
146 : * end. For nicely shaped elements, the containing Elem is often the
147 : * first or second one found in the Nanoflann results.
148 : *
149 : * 2.) I'm not sure about the relative cost of requesting more
150 : * results from the Nanoflann search than one actually needs, but
151 : * presumably reducing _num_results will result in better
152 : * performance for your particular application up to a point. If, on
153 : * the other hand, _num_results is too small, then you risk
154 : * having a false negative result, so take this into account when
155 : * choosing the parameter for a particular application.
156 : */
157 : std::size_t get_num_results() const;
158 : void set_num_results(std::size_t val);
159 :
160 : //
161 : // Required Nanoflann typedefs and APIs
162 : //
163 :
164 : /**
165 : * Floating point type used for storing coordinates
166 : */
167 : typedef Real coord_t;
168 :
169 : /**
170 : * Must return the number of data points
171 : */
172 : std::size_t kdtree_get_point_count() const;
173 :
174 : /**
175 : * Returns the squared distance between the vector (p1[0], p1[1], p1[2])
176 : * and the vertex average of Elem "idx_p2" stored in _mesh
177 : */
178 : coord_t kdtree_distance(const coord_t * p1,
179 : const std::size_t idx_p2,
180 : std::size_t size) const;
181 :
182 : /**
183 : * Returns the dim'th component of the idx'th vertex average.
184 : */
185 : coord_t kdtree_get_pt(const std::size_t idx, int dim) const;
186 :
187 : /**
188 : * Optional bounding-box computation: return false to default to a
189 : * standard bbox computation loop. Return true if the BBOX can
190 : * be computed more efficiently (and returned in "bb") than the
191 : * standard bounding box computation. The BBOX template parameter
192 : * must at least support bb.size() to find out the expected
193 : * dimensionality of the box (e.g. 2 or 3 for point clouds).
194 : */
195 : template <class BBOX>
196 0 : bool kdtree_get_bbox(BBOX & /*bb*/) const { return false; }
197 :
198 : protected:
199 :
200 : /**
201 : * \p true if out-of-mesh mode is enabled. See \p
202 : * enable_out_of_mesh_mode() for details.
203 : */
204 : bool _out_of_mesh_mode;
205 :
206 : /**
207 : * The number of results returned by Nanoflann when operator() is
208 : * called.
209 : */
210 : std::size_t _num_results;
211 :
212 : /**
213 : * Lists of Points and ids which make up the "point cloud" that is
214 : * to be searched via Nanoflann. The point cloud can be comprised of
215 : * Elem vertex averages or mesh Nodes, depending on the type of tree
216 : * created. We keep two separate vectors since the Points in the
217 : * cloud may not be numbered contiguously in general.
218 : *
219 : * These are shared_ptrs to vectors since, if we are not the "master"
220 : * PointLocator, they need to point at the master's vectors instead.
221 : */
222 : std::shared_ptr<std::vector<const Elem *>> _elems;
223 : std::shared_ptr<std::vector<Point>> _point_cloud;
224 :
225 : // kd_tree will be initialized during init() and then automatically
226 : // cleaned up by the destructor. We always create a LIBMESH_DIM
227 : // dimensional Nanoflann object.
228 : typedef nanoflann::L2_Simple_Adaptor<Real, PointLocatorNanoflann> adapter_t;
229 : typedef nanoflann::KDTreeSingleIndexAdaptor<adapter_t, PointLocatorNanoflann, LIBMESH_DIM> kd_tree_t;
230 : std::shared_ptr<kd_tree_t> _kd_tree;
231 :
232 : /**
233 : * The operator() functions on PointLocator-derived classes are
234 : * const, so to avoid re-allocating these result data structures every
235 : * time operator() is called, they have to be mutable.
236 : */
237 : mutable std::vector<std::size_t> _ret_index;
238 : mutable std::vector<Real> _out_dist_sqr;
239 :
240 : /**
241 : * Vector of indices used by indirect sort. Stored as a class member
242 : * so we don't have to allocate it from scratch for every search.
243 : */
244 : mutable std::vector<std::size_t> _b;
245 :
246 : /**
247 : * Helper function that wraps the call to the KDTree's
248 : * findNeighbors() routine. Must be passed the Point to search for
249 : * and the number of results to return. Stores the results of the
250 : * search in the _ret_index and _out_dist_sqr class members and
251 : * returns a KNNResultSet, which has pointers to the index and
252 : * distance data.
253 : */
254 : nanoflann::KNNResultSet<Real>
255 : kd_tree_find_neighbors(const Point & p,
256 : std::size_t num_results) const;
257 : };
258 :
259 : } // namespace libMesh
260 :
261 : #endif // LIBMESH_HAVE_NANOFLANN
262 : #endif // LIBMESH_POINT_LOCATOR_NANOFLANN_H
|
