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 : // C++ includes
21 :
22 : // Local Includes
23 : #include "libmesh/elem.h"
24 : #include "libmesh/libmesh_logging.h"
25 : #include "libmesh/mesh_base.h"
26 : #include "libmesh/mesh_tools.h"
27 : #include "libmesh/point_locator_tree.h"
28 : #include "libmesh/tree.h"
29 :
30 : namespace libMesh
31 : {
32 :
33 :
34 :
35 : //------------------------------------------------------------------
36 : // PointLocator methods
37 6 : PointLocatorTree::PointLocatorTree (const MeshBase & mesh,
38 6 : const PointLocatorBase * master) :
39 : PointLocatorBase (mesh,master),
40 : _tree (nullptr),
41 2 : _element (nullptr),
42 2 : _out_of_mesh_mode(false),
43 2 : _target_bin_size (200),
44 6 : _build_type(Trees::NODES)
45 : {
46 6 : this->init(_build_type);
47 6 : }
48 :
49 :
50 :
51 500073 : PointLocatorTree::PointLocatorTree (const MeshBase & mesh,
52 : const Trees::BuildType build_type,
53 500073 : const PointLocatorBase * master) :
54 : PointLocatorBase (mesh,master),
55 : _tree (nullptr),
56 387479 : _element (nullptr),
57 387479 : _out_of_mesh_mode(false),
58 387479 : _target_bin_size (200),
59 500073 : _build_type(build_type)
60 : {
61 500073 : this->init(_build_type);
62 500073 : }
63 :
64 :
65 :
66 1162441 : PointLocatorTree::~PointLocatorTree () = default;
67 :
68 :
69 :
70 0 : void PointLocatorTree::clear ()
71 : {
72 0 : this->_initialized = false;
73 :
74 : // reset() actually frees the memory if we are master, otherwise it
75 : // just reduces the ref. count.
76 0 : _tree.reset();
77 0 : }
78 :
79 :
80 :
81 0 : void PointLocatorTree::init()
82 : {
83 0 : this->init(_build_type);
84 0 : }
85 :
86 :
87 :
88 500079 : void PointLocatorTree::init (Trees::BuildType build_type)
89 : {
90 56499 : libmesh_assert (!this->_tree);
91 :
92 500079 : if (this->_initialized)
93 : {
94 : // Warn that we are already initialized
95 0 : libMesh::err << "Warning: PointLocatorTree already initialized! Will ignore this call..." << std::endl;
96 :
97 : // Further warn if we try to init() again with a different build_type
98 0 : if (_build_type != build_type)
99 : {
100 0 : libMesh::err << "Warning: PointLocatorTree is using build_type = " << _build_type << ".\n"
101 0 : << "Your requested build_type, " << build_type << " will not be used!" << std::endl;
102 : }
103 : }
104 :
105 : else
106 : {
107 : // Let the requested build_type override the _build_type we were
108 : // constructed with. This is no big deal since we have not been
109 : // initialized before.
110 500079 : _build_type = build_type;
111 :
112 500079 : if (this->_master == nullptr)
113 : {
114 1012 : LOG_SCOPE("init(no master)", "PointLocatorTree");
115 :
116 : if (LIBMESH_DIM == 1)
117 : _tree = std::make_shared<Trees::BinaryTree>(this->_mesh, get_target_bin_size(), _build_type);
118 : else if (LIBMESH_DIM == 2)
119 : _tree = std::make_shared<Trees::QuadTree>(this->_mesh, get_target_bin_size(), _build_type);
120 12547 : else if (this->_mesh.mesh_dimension() == 3) // && LIBMESH_DIM==3
121 5884 : _tree = std::make_shared<Trees::OctTree>(this->_mesh, get_target_bin_size(), _build_type);
122 : else
123 : {
124 : // LIBMESH_DIM==3 but we have a mesh with only 1D/2D
125 : // elements, which needs special consideration. If the
126 : // mesh is planar XY, we want to build a QuadTree to
127 : // search efficiently. If the mesh is truly a manifold,
128 : // then we need an octree
129 422 : bool is_planar_xy = false;
130 :
131 : // Build the bounding box for the mesh. If the delta-z bound is
132 : // negligibly small then we can use a quadtree.
133 9563 : BoundingBox bbox = MeshTools::create_bounding_box(this->_mesh);
134 :
135 : const Real
136 9563 : Dx = bbox.second(0) - bbox.first(0),
137 9563 : Dz = bbox.second(2) - bbox.first(2);
138 :
139 : // In order to satisfy is_planar_xy the mesh should be planar and should
140 : // also be in the z=0 plane, since otherwise it is incorrect to use a
141 : // QuadTree since QuadTrees assume z=0.
142 9985 : if ( (std::abs(Dz/(Dx + 1.e-20)) < 1e-10) && (std::abs(bbox.second(2)) < 1.e-10) )
143 418 : is_planar_xy = true;
144 :
145 422 : if (is_planar_xy)
146 18424 : _tree = std::make_shared<Trees::QuadTree>(this->_mesh, get_target_bin_size(), _build_type);
147 : else
148 280 : _tree = std::make_shared<Trees::OctTree>(this->_mesh, get_target_bin_size(), _build_type);
149 : }
150 : }
151 :
152 : else
153 : {
154 : // We are _not_ the master. Let our Tree point to
155 : // the master's tree. But for this we first transform
156 : // the master in a state for which we are friends.
157 : // And make sure the master has a tree!
158 : const PointLocatorTree * my_master =
159 55993 : cast_ptr<const PointLocatorTree *>(this->_master);
160 :
161 487532 : if (my_master->initialized())
162 55993 : this->_tree = my_master->_tree;
163 : else
164 0 : libmesh_error_msg("ERROR: Initialize master first, then servants!");
165 : }
166 :
167 : // Not all PointLocators may own a tree, but all of them
168 : // use their own element pointer. Let the element pointer
169 : // be unique for every interpolator.
170 : // Suppose the interpolators are used concurrently
171 : // at different locations in the mesh, then it makes quite
172 : // sense to have unique start elements.
173 500079 : this->_element = nullptr;
174 : }
175 :
176 : // ready for take-off
177 500079 : this->_initialized = true;
178 500079 : }
179 :
180 4532083 : const Elem * PointLocatorTree::operator() (const Point & p,
181 : const std::set<subdomain_id_type> * allowed_subdomains) const
182 : {
183 327139 : libmesh_assert (this->_initialized);
184 :
185 654278 : LOG_SCOPE("operator()", "PointLocatorTree");
186 :
187 : // If we're provided with an allowed_subdomains list and have a cached element, make sure it complies
188 7365686 : if (allowed_subdomains && this->_element && !allowed_subdomains->count(this->_element->subdomain_id()))
189 1511 : this->_element = nullptr;
190 :
191 4532083 : if (this->_element != nullptr)
192 : {
193 : // If the user specified a custom tolerance, we actually call
194 : // Elem::close_to_point() instead, since Elem::contains_point()
195 : // warns about using non-default BoundingBox tolerances.
196 4213450 : if (_use_contains_point_tol && !(this->_element->close_to_point(p, _contains_point_tol)))
197 62172 : this->_element = nullptr;
198 :
199 : // Otherwise, just call contains_point(p) with default tolerances.
200 4151278 : else if (!(this->_element->contains_point(p)))
201 337686 : this->_element = nullptr;
202 : }
203 :
204 : // First check the element from last time before asking the tree
205 4532083 : if (this->_element==nullptr)
206 : {
207 : // ask the tree
208 718491 : if (_use_contains_point_tol)
209 62551 : this->_element = this->_tree->find_element(p, allowed_subdomains, _contains_point_tol);
210 : else
211 655940 : this->_element = this->_tree->find_element(p, allowed_subdomains);
212 :
213 718491 : if (this->_element == nullptr)
214 : {
215 : // If we haven't found the element, we may want to do a linear
216 : // search using a tolerance.
217 152101 : if (_use_close_to_point_tol)
218 : {
219 0 : if (_verbose)
220 : {
221 0 : libMesh::out << "Performing linear search using close-to-point tolerance "
222 0 : << _close_to_point_tol
223 0 : << std::endl;
224 : }
225 :
226 0 : this->_element =
227 0 : this->perform_linear_search(p,
228 : allowed_subdomains,
229 : /*use_close_to_point*/ true,
230 0 : _close_to_point_tol);
231 :
232 0 : return this->_element;
233 : }
234 :
235 : // No element seems to contain this point. In theory, our
236 : // tree now correctly handles curved elements. In
237 : // out-of-mesh mode this is sometimes expected, and we can
238 : // just return nullptr without searching further. Out of
239 : // out-of-mesh mode, something must have gone wrong.
240 83 : libmesh_assert_equal_to (_out_of_mesh_mode, true);
241 :
242 83 : return this->_element;
243 : }
244 : }
245 :
246 : // If we found an element, it should be active
247 327056 : libmesh_assert (!this->_element || this->_element->active());
248 :
249 : // If we found an element and have a restriction list, they better match
250 327056 : libmesh_assert (!this->_element || !allowed_subdomains || allowed_subdomains->count(this->_element->subdomain_id()));
251 :
252 : // return the element
253 4379982 : return this->_element;
254 : }
255 :
256 :
257 256693 : void PointLocatorTree::operator() (const Point & p,
258 : std::set<const Elem *> & candidate_elements,
259 : const std::set<subdomain_id_type> * allowed_subdomains) const
260 : {
261 32252 : libmesh_assert (this->_initialized);
262 :
263 64504 : LOG_SCOPE("operator() - Version 2", "PointLocatorTree");
264 :
265 : // ask the tree
266 256693 : this->_tree->find_elements (p, candidate_elements, allowed_subdomains, _close_to_point_tol);
267 256693 : }
268 :
269 :
270 :
271 0 : const Elem * PointLocatorTree::perform_linear_search(const Point & p,
272 : const std::set<subdomain_id_type> * allowed_subdomains,
273 : bool use_close_to_point,
274 : Real close_to_point_tolerance) const
275 : {
276 0 : LOG_SCOPE("perform_linear_search", "PointLocatorTree");
277 :
278 : // The type of iterator depends on the Trees::BuildType
279 : // used for this PointLocator. If it's
280 : // TREE_LOCAL_ELEMENTS, we only want to double check
281 : // local elements during this linear search.
282 : SimpleRange<MeshBase::const_element_iterator> r =
283 0 : this->_build_type == Trees::LOCAL_ELEMENTS ?
284 0 : this->_mesh.active_local_element_ptr_range() :
285 0 : this->_mesh.active_element_ptr_range();
286 :
287 0 : for (const auto & elem : r)
288 : {
289 0 : if (!allowed_subdomains ||
290 0 : allowed_subdomains->count(elem->subdomain_id()))
291 : {
292 0 : if (!use_close_to_point)
293 : {
294 0 : if (elem->contains_point(p, _contains_point_tol))
295 0 : return elem;
296 : }
297 : else
298 : {
299 0 : if (elem->close_to_point(p, close_to_point_tolerance))
300 0 : return elem;
301 : }
302 : }
303 : }
304 :
305 0 : return nullptr;
306 0 : }
307 :
308 :
309 0 : std::set<const Elem *> PointLocatorTree::perform_fuzzy_linear_search(const Point & p,
310 : const std::set<subdomain_id_type> * allowed_subdomains,
311 : Real close_to_point_tolerance) const
312 : {
313 0 : LOG_SCOPE("perform_fuzzy_linear_search", "PointLocatorTree");
314 :
315 0 : std::set<const Elem *> candidate_elements;
316 :
317 : // The type of iterator depends on the Trees::BuildType
318 : // used for this PointLocator. If it's
319 : // TREE_LOCAL_ELEMENTS, we only want to double check
320 : // local elements during this linear search.
321 : SimpleRange<MeshBase::const_element_iterator> r =
322 0 : this->_build_type == Trees::LOCAL_ELEMENTS ?
323 0 : this->_mesh.active_local_element_ptr_range() :
324 0 : this->_mesh.active_element_ptr_range();
325 :
326 0 : for (const auto & elem : r)
327 0 : if ((!allowed_subdomains || allowed_subdomains->count(elem->subdomain_id())) && elem->close_to_point(p, close_to_point_tolerance))
328 0 : candidate_elements.insert(elem);
329 :
330 0 : return candidate_elements;
331 0 : }
332 :
333 :
334 :
335 246870 : void PointLocatorTree::enable_out_of_mesh_mode ()
336 : {
337 : // Out-of-mesh mode should now work properly even on meshes with
338 : // non-affine elements.
339 246870 : _out_of_mesh_mode = true;
340 246870 : }
341 :
342 :
343 0 : void PointLocatorTree::disable_out_of_mesh_mode ()
344 : {
345 0 : _out_of_mesh_mode = false;
346 0 : }
347 :
348 :
349 0 : void PointLocatorTree::set_target_bin_size (unsigned int target_bin_size)
350 : {
351 0 : _target_bin_size = target_bin_size;
352 0 : }
353 :
354 :
355 12547 : unsigned int PointLocatorTree::get_target_bin_size () const
356 : {
357 12547 : return _target_bin_size;
358 : }
359 :
360 :
361 : } // namespace libMesh
|
