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 : #include <set>
22 : #include <array>
23 :
24 : // Local includes
25 : #include "libmesh/libmesh_config.h"
26 : #include "libmesh/tree_node.h"
27 : #include "libmesh/mesh_base.h"
28 : #include "libmesh/elem.h"
29 :
30 : namespace libMesh
31 : {
32 :
33 : // ------------------------------------------------------------
34 : // TreeNode class methods
35 : template <unsigned int N>
36 126430 : bool TreeNode<N>::insert (const Node * nd)
37 : {
38 49762 : libmesh_assert(nd);
39 49762 : libmesh_assert_less (nd->id(), mesh.n_nodes());
40 :
41 : // Return if we don't bound the node
42 126430 : if (!this->bounds_node(nd))
43 38928 : return false;
44 :
45 : // Add the node to ourself if we are active
46 27583 : if (this->active())
47 : {
48 21813 : nodes.push_back (nd);
49 :
50 : // Refine ourself if we reach the target bin size for a TreeNode.
51 30349 : if (nodes.size() == tgt_bin_size)
52 46 : this->refine();
53 :
54 21813 : return true;
55 : }
56 :
57 : // If we are not active simply pass the node along to
58 : // our children
59 2298 : libmesh_assert_equal_to (children.size(), N);
60 :
61 2298 : bool was_inserted = false;
62 51930 : for (unsigned int c=0; c<N; c++)
63 64544 : if (children[c]->insert (nd))
64 2962 : was_inserted = true;
65 2298 : return was_inserted;
66 : }
67 :
68 :
69 :
70 : template <unsigned int N>
71 21050651 : bool TreeNode<N>::insert (const Elem * elem)
72 : {
73 2184147 : libmesh_assert(elem);
74 :
75 : // We first want to find the corners of the cuboid surrounding the cell.
76 21050651 : const BoundingBox bbox = elem->loose_bounding_box();
77 :
78 : // If we are using a QuadTree, it's either because LIBMESH_DIM==2 or
79 : // we have a planar xy mesh. Either way, the bounding box
80 : // comparison in this case needs to do something slightly different
81 : // for the z-coordinate.
82 2184147 : bool bboxes_intersect = false;
83 :
84 : if (N == 8) // OctTree
85 7682413 : bboxes_intersect = this->bounding_box.intersects(bbox);
86 : else if (N == 4) // QuadTree
87 : {
88 : // Perform a specialized BoundingBox intersection check that
89 : // ignores z-coords. Copied from geom/bounding_box.C Check for
90 : // "real" intersection in the x and y directions, then check
91 : // that the z-coordinate is "close".
92 :
93 : // Helper macro
94 : #define IS_BETWEEN(min, check, max) \
95 : ((min) <= (check) && (check) <= (max))
96 :
97 : // Make local variables first to make things more clear in a moment
98 1954713 : const Real & elem_min_x = bbox.first(0);
99 1954713 : const Real & elem_max_x = bbox.second(0);
100 1954713 : const Real & tree_min_x = this->bounding_box.first(0);
101 1954713 : const Real & tree_max_x = this->bounding_box.second(0);
102 :
103 1954713 : const Real & elem_min_y = bbox.first(1);
104 1954713 : const Real & elem_max_y = bbox.second(1);
105 1954713 : const Real & tree_min_y = this->bounding_box.first(1);
106 1954713 : const Real & tree_max_y = this->bounding_box.second(1);
107 :
108 1954713 : bool x_int =
109 3420557 : IS_BETWEEN(elem_min_x, tree_min_x, elem_max_x) ||
110 12775415 : IS_BETWEEN(elem_min_x, tree_max_x, elem_max_x) ||
111 24446093 : IS_BETWEEN(tree_min_x, elem_min_x, tree_max_x) ||
112 5458516 : IS_BETWEEN(tree_min_x, elem_max_x, tree_max_x);
113 :
114 1954713 : bool y_int =
115 9723888 : IS_BETWEEN(elem_min_y, tree_min_y, elem_max_y) ||
116 4464898 : IS_BETWEEN(elem_min_y, tree_max_y, elem_max_y) ||
117 18181875 : IS_BETWEEN(tree_min_y, elem_min_y, tree_max_y) ||
118 1506212 : IS_BETWEEN(tree_min_y, elem_max_y, tree_max_y);
119 :
120 : // When LIBMESH_DIM==3, check that the z-coordinates of the elem
121 : // bbox and the tree bbox are "close" since the QuadTree is
122 : // meant to work in the case when the mesh is planar_xy but not
123 : // necessarily lying in the z=0 plane.
124 1954713 : bool z_match = true;
125 : if (LIBMESH_DIM == 3)
126 : {
127 1954713 : const Real & elem_min_z = bbox.first(2);
128 1954713 : const Real & elem_max_z = bbox.second(2);
129 1954713 : const Real & tree_min_z = this->bounding_box.first(2);
130 1954713 : const Real & tree_max_z = this->bounding_box.second(2);
131 :
132 1954713 : z_match =
133 15093517 : (std::abs(elem_min_z - elem_max_z) < TOLERANCE) &&
134 15093517 : (std::abs(tree_min_z - tree_max_z) < TOLERANCE) &&
135 13138804 : (std::abs(elem_min_z - tree_max_z) < TOLERANCE);
136 : }
137 :
138 13138804 : bboxes_intersect = z_match && x_int && y_int;
139 : }
140 : else // binary tree
141 : {
142 : // TODO: implement 1D bounding box intersection check
143 0 : libmesh_not_implemented();
144 : }
145 :
146 : // Next, find out whether this cuboid has got non-empty intersection
147 : // with the bounding box of the current tree node.
148 : //
149 : // If not, we should not care about this element.
150 9866560 : if (!bboxes_intersect)
151 1375064 : return false;
152 :
153 : // Only add the element if we are active
154 8345119 : if (this->active())
155 : {
156 6385051 : elements.push_back (elem);
157 :
158 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
159 :
160 : // flag indicating this node contains
161 : // infinite elements
162 1376368 : if (elem->infinite())
163 0 : this->contains_ifems = true;
164 :
165 : #endif
166 :
167 1142814 : unsigned int element_count = cast_int<unsigned int>(elements.size());
168 6385051 : if (!mesh.get_count_lower_dim_elems_in_point_locator())
169 : {
170 0 : const std::set<unsigned char> & elem_dimensions = mesh.elem_dimensions();
171 0 : if (elem_dimensions.size() > 1)
172 : {
173 0 : element_count = 0;
174 0 : unsigned char highest_dim_elem = *elem_dimensions.rbegin();
175 0 : for (const Elem * other_elem : elements)
176 0 : if (other_elem->dim() == highest_dim_elem)
177 0 : element_count++;
178 : }
179 : }
180 :
181 : // Refine ourself if we reach the target bin size for a TreeNode.
182 6385051 : if (element_count == tgt_bin_size)
183 10120 : this->refine();
184 :
185 6385051 : return true;
186 : }
187 :
188 : // If we are not active simply pass the element along to
189 : // our children
190 237676 : libmesh_assert_equal_to (children.size(), N);
191 :
192 237676 : bool was_inserted = false;
193 12158648 : for (unsigned int c=0; c<N; c++)
194 11217732 : if (children[c]->insert (elem))
195 306008 : was_inserted = true;
196 237676 : return was_inserted;
197 : }
198 :
199 :
200 :
201 : template <unsigned int N>
202 10166 : void TreeNode<N>::refine ()
203 : {
204 : // Huh? better be active...
205 1198 : libmesh_assert (this->active());
206 1198 : libmesh_assert (children.empty());
207 :
208 : // A TreeNode<N> has by definition N children
209 10166 : children.resize(N);
210 :
211 : // Scale up the target bin size in child TreeNodes if we have reached
212 : // the maximum number of refinement levels.
213 10166 : unsigned int new_target_bin_size = tgt_bin_size;
214 10166 : if (level() >= target_bin_size_increase_level)
215 : {
216 0 : new_target_bin_size *= 2;
217 : }
218 :
219 57702 : for (unsigned int c=0; c<N; c++)
220 : {
221 : // Create the child and set its bounding box.
222 47536 : children[c] = std::make_unique<TreeNode<N>>(mesh, new_target_bin_size, this);
223 52592 : children[c]->set_bounding_box(this->create_bounding_box(c));
224 :
225 : // Pass off our nodes to our children
226 121136 : for (const Node * node : nodes)
227 102400 : children[c]->insert(node);
228 :
229 : // Pass off our elements to our children
230 9481136 : for (const Elem * elem : elements)
231 10416000 : children[c]->insert(elem);
232 : }
233 :
234 : // We don't need to store nodes or elements any more, they have been
235 : // added to the children. Use the "swap trick" to actually reduce
236 : // the capacity of these vectors.
237 8968 : std::vector<const Node *>().swap(nodes);
238 8968 : std::vector<const Elem *>().swap(elements);
239 :
240 1198 : libmesh_assert_equal_to (nodes.capacity(), 0);
241 1198 : libmesh_assert_equal_to (elements.capacity(), 0);
242 10166 : }
243 :
244 :
245 :
246 : template <unsigned int N>
247 60083 : void TreeNode<N>::set_bounding_box (const std::pair<Point, Point> & bbox)
248 : {
249 60083 : bounding_box = bbox;
250 60083 : }
251 :
252 :
253 :
254 : template <unsigned int N>
255 99524 : bool TreeNode<N>::bounds_node (const Node * nd,
256 : Real relative_tol) const
257 : {
258 49762 : libmesh_assert(nd);
259 126430 : return bounds_point(*nd, relative_tol);
260 : }
261 :
262 :
263 :
264 : template <unsigned int N>
265 2003097 : bool TreeNode<N>::bounds_point (const Point & p,
266 : Real relative_tol) const
267 : {
268 320866 : const Point & min = bounding_box.first;
269 320866 : const Point & max = bounding_box.second;
270 :
271 2003097 : const Real tol = (max - min).norm() * relative_tol;
272 :
273 2003097 : if ((p(0) >= min(0) - tol)
274 1785889 : && (p(0) <= max(0) + tol)
275 : #if LIBMESH_DIM > 1
276 1523895 : && (p(1) >= min(1) - tol)
277 1425245 : && (p(1) <= max(1) + tol)
278 : #endif
279 : #if LIBMESH_DIM > 2
280 1282983 : && (p(2) >= min(2) - tol)
281 3401820 : && (p(2) <= max(2) + tol)
282 : #endif
283 : )
284 1209691 : return true;
285 :
286 201254 : return false;
287 : }
288 :
289 :
290 :
291 : template <unsigned int N>
292 : BoundingBox
293 47536 : TreeNode<N>::create_bounding_box (unsigned int c) const
294 : {
295 : switch (N)
296 : {
297 : // How to refine an OctTree Node
298 : case 8:
299 : {
300 13744 : const Real xmin = bounding_box.first(0);
301 13744 : const Real ymin = bounding_box.first(1);
302 13744 : const Real zmin = bounding_box.first(2);
303 :
304 13744 : const Real xmax = bounding_box.second(0);
305 13744 : const Real ymax = bounding_box.second(1);
306 13744 : const Real zmax = bounding_box.second(2);
307 :
308 13744 : const Real xc = .5*(xmin + xmax);
309 13744 : const Real yc = .5*(ymin + ymax);
310 13744 : const Real zc = .5*(zmin + zmax);
311 :
312 13744 : switch (c)
313 : {
314 66 : case 0:
315 : return BoundingBox (Point(xmin, ymin, zmin),
316 66 : Point(xc, yc, zc));
317 66 : case 1:
318 : return BoundingBox (Point(xc, ymin, zmin),
319 66 : Point(xmax, yc, zc));
320 66 : case 2:
321 : return BoundingBox (Point(xmin, yc, zmin),
322 66 : Point(xc, ymax, zc));
323 66 : case 3:
324 : return BoundingBox (Point(xc, yc, zmin),
325 66 : Point(xmax, ymax, zc));
326 66 : case 4:
327 : return BoundingBox (Point(xmin, ymin, zc),
328 66 : Point(xc, yc, zmax));
329 66 : case 5:
330 : return BoundingBox (Point(xc, ymin, zc),
331 66 : Point(xmax, yc, zmax));
332 66 : case 6:
333 : return BoundingBox (Point(xmin, yc, zc),
334 66 : Point(xc, ymax, zmax));
335 66 : case 7:
336 : return BoundingBox (Point(xc, yc, zc),
337 66 : Point(xmax, ymax, zmax));
338 0 : default:
339 0 : libmesh_error_msg("c >= N! : " << c);
340 : }
341 :
342 : break;
343 : } // case 8
344 :
345 : // How to refine an QuadTree Node
346 : case 4:
347 : {
348 33792 : const Real xmin = bounding_box.first(0);
349 33792 : const Real ymin = bounding_box.first(1);
350 :
351 33792 : const Real xmax = bounding_box.second(0);
352 33792 : const Real ymax = bounding_box.second(1);
353 :
354 33792 : const Real xc = .5*(xmin + xmax);
355 33792 : const Real yc = .5*(ymin + ymax);
356 :
357 33792 : switch (c)
358 : {
359 1132 : case 0:
360 : return BoundingBox (Point(xmin, ymin),
361 1132 : Point(xc, yc));
362 1132 : case 1:
363 : return BoundingBox (Point(xc, ymin),
364 1132 : Point(xmax, yc));
365 1132 : case 2:
366 : return BoundingBox (Point(xmin, yc),
367 1132 : Point(xc, ymax));
368 1132 : case 3:
369 : return BoundingBox (Point(xc, yc),
370 1132 : Point(xmax, ymax));
371 0 : default:
372 0 : libmesh_error_msg("c >= N!");
373 : }
374 :
375 : break;
376 : } // case 4
377 :
378 : // How to refine a BinaryTree Node
379 : case 2:
380 : {
381 0 : const Real xmin = bounding_box.first(0);
382 :
383 0 : const Real xmax = bounding_box.second(0);
384 :
385 0 : const Real xc = .5*(xmin + xmax);
386 :
387 0 : switch (c)
388 : {
389 0 : case 0:
390 : return BoundingBox (Point(xmin),
391 0 : Point(xc));
392 0 : case 1:
393 : return BoundingBox (Point(xc),
394 0 : Point(xmax));
395 0 : default:
396 0 : libmesh_error_msg("c >= N!");
397 : }
398 :
399 : break;
400 : } // case 2
401 :
402 : default:
403 : libmesh_error_msg("Only implemented for Octrees, QuadTrees, and Binary Trees!");
404 : }
405 : }
406 :
407 :
408 :
409 : template <unsigned int N>
410 0 : void TreeNode<N>::print_nodes(std::ostream & out_stream) const
411 : {
412 0 : if (this->active())
413 : {
414 0 : out_stream << "TreeNode Level: " << this->level() << std::endl;
415 :
416 0 : for (const Node * node : nodes)
417 0 : out_stream << " " << node->id();
418 :
419 0 : out_stream << std::endl << std::endl;
420 : }
421 : else
422 0 : for (const auto & child : children)
423 0 : child->print_nodes();
424 0 : }
425 :
426 :
427 :
428 : template <unsigned int N>
429 0 : void TreeNode<N>::print_elements(std::ostream & out_stream) const
430 : {
431 0 : if (this->active())
432 : {
433 0 : out_stream << "TreeNode Level: " << this->level() << std::endl;
434 :
435 0 : for (const auto & elem : elements)
436 0 : out_stream << " " << elem;
437 :
438 0 : out_stream << std::endl << std::endl;
439 : }
440 : else
441 0 : for (const auto & child : children)
442 0 : child->print_elements();
443 0 : }
444 :
445 :
446 :
447 : template <unsigned int N>
448 0 : void TreeNode<N>::transform_nodes_to_elements (std::vector<std::vector<const Elem *>> & nodes_to_elem)
449 : {
450 0 : if (this->active())
451 : {
452 0 : elements.clear();
453 :
454 : // Temporarily use a set. Since multiple nodes
455 : // will likely map to the same element we use a
456 : // set to eliminate the duplication.
457 0 : std::set<const Elem *> elements_set;
458 :
459 0 : for (const Node * node : nodes)
460 : {
461 : // the actual global node number we are replacing
462 : // with the connected elements
463 0 : const dof_id_type node_number = node->id();
464 :
465 0 : libmesh_assert_less (node_number, mesh.n_nodes());
466 0 : libmesh_assert_less (node_number, nodes_to_elem.size());
467 :
468 0 : for (const Elem * elem : nodes_to_elem[node_number])
469 0 : elements_set.insert(elem);
470 : }
471 :
472 : // Done with the nodes.
473 0 : std::vector<const Node *>().swap(nodes);
474 :
475 : // Now the set is built. We can copy this to the
476 : // vector. Note that the resulting vector will
477 : // already be sorted, and will require less memory
478 : // than the set.
479 0 : elements.reserve(elements_set.size());
480 :
481 0 : for (const auto & elem : elements_set)
482 : {
483 0 : elements.push_back(elem);
484 :
485 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
486 :
487 : // flag indicating this node contains
488 : // infinite elements
489 0 : if (elem->infinite())
490 0 : this->contains_ifems = true;
491 :
492 : #endif
493 : }
494 : }
495 : else
496 0 : for (auto & child : children)
497 0 : child->transform_nodes_to_elements (nodes_to_elem);
498 0 : }
499 :
500 :
501 :
502 : template <unsigned int N>
503 374 : void TreeNode<N>::transform_nodes_to_elements (std::unordered_map<dof_id_type, std::vector<const Elem *>> & nodes_to_elem)
504 : {
505 374 : if (this->active())
506 : {
507 328 : elements.clear();
508 :
509 : // Temporarily use a set. Since multiple nodes
510 : // will likely map to the same element we use a
511 : // set to eliminate the duplication.
512 256 : std::set<const Elem *> elements_set;
513 :
514 12941 : for (const Node * node : nodes)
515 : {
516 : // the actual global node number we are replacing
517 : // with the connected elements
518 12613 : const dof_id_type node_number = node->id();
519 :
520 4936 : libmesh_assert_less (node_number, mesh.n_nodes());
521 :
522 4936 : auto & my_elems = nodes_to_elem[node_number];
523 7677 : elements_set.insert(my_elems.begin(), my_elems.end());
524 : }
525 :
526 : // Done with the nodes.
527 200 : std::vector<const Node *>().swap(nodes);
528 :
529 : // Now the set is built. We can copy this to the
530 : // vector. Note that the resulting vector will
531 : // already be sorted, and will require less memory
532 : // than the set.
533 328 : elements.reserve(elements_set.size());
534 :
535 9238 : for (const auto & elem : elements_set)
536 : {
537 8910 : elements.push_back(elem);
538 :
539 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
540 :
541 : // flag indicating this node contains
542 : // infinite elements
543 8910 : if (elem->infinite())
544 5197 : this->contains_ifems = true;
545 :
546 : #endif
547 : }
548 : }
549 : else
550 414 : for (auto & child : children)
551 368 : child->transform_nodes_to_elements (nodes_to_elem);
552 374 : }
553 :
554 :
555 :
556 : template <unsigned int N>
557 0 : unsigned int TreeNode<N>::n_active_bins() const
558 : {
559 0 : if (this->active())
560 0 : return 1;
561 :
562 : else
563 : {
564 0 : unsigned int sum=0;
565 :
566 0 : for (const auto & child : children)
567 0 : sum += child->n_active_bins();
568 :
569 0 : return sum;
570 : }
571 : }
572 :
573 :
574 :
575 : template <unsigned int N>
576 : const Elem *
577 810649 : TreeNode<N>::find_element (const Point & p,
578 : const std::set<subdomain_id_type> * allowed_subdomains,
579 : Real relative_tol) const
580 : {
581 810649 : if (this->active())
582 : {
583 : // Only check our children if the point is in our bounding box
584 : // or if the node contains infinite elements
585 718521 : if (this->bounds_point(p, relative_tol) || this->contains_ifems)
586 : // Search the active elements in the active TreeNode.
587 25951418 : for (const auto & elem : elements)
588 25799317 : if (!allowed_subdomains || allowed_subdomains->count(elem->subdomain_id()))
589 25078350 : if (elem->active())
590 : {
591 1826130 : bool found = relative_tol > TOLERANCE
592 23270970 : ? elem->close_to_point(p, relative_tol)
593 24956755 : : elem->contains_point(p);
594 :
595 25078236 : if (found)
596 566390 : return elem;
597 : }
598 :
599 : // The point was not found in any element
600 152131 : return nullptr;
601 : }
602 : else
603 77012 : return this->find_element_in_children(p,allowed_subdomains,
604 92128 : relative_tol);
605 : }
606 :
607 :
608 :
609 : template <unsigned int N>
610 : void
611 368049 : TreeNode<N>::find_elements (const Point & p,
612 : std::set<const Elem *> & candidate_elements,
613 : const std::set<subdomain_id_type> * allowed_subdomains,
614 : Real relative_tol) const
615 : {
616 368049 : if (this->active())
617 : {
618 : // Only check our children if the point is in our bounding box
619 : // or if the node contains infinite elements
620 260133 : if (this->bounds_point(p, relative_tol) || this->contains_ifems)
621 : // Search the active elements in the active TreeNode.
622 15054371 : for (const auto & elem : elements)
623 14794238 : if (!allowed_subdomains || allowed_subdomains->count(elem->subdomain_id()))
624 13985846 : if (elem->active())
625 : {
626 4301842 : bool found = relative_tol > TOLERANCE
627 9684004 : ? elem->close_to_point(p, relative_tol)
628 13985846 : : elem->contains_point(p);
629 :
630 13985846 : if (found)
631 465126 : candidate_elements.insert(elem);
632 : }
633 : }
634 : else
635 107916 : this->find_elements_in_children(p, candidate_elements,
636 : allowed_subdomains, relative_tol);
637 368049 : }
638 :
639 :
640 :
641 : template <unsigned int N>
642 92128 : const Elem * TreeNode<N>::find_element_in_children (const Point & p,
643 : const std::set<subdomain_id_type> * allowed_subdomains,
644 : Real relative_tol) const
645 : {
646 7558 : libmesh_assert (!this->active());
647 :
648 : // value-initialization sets all array members to false
649 99022 : auto searched_child = std::array<bool, N>();
650 :
651 : // First only look in the children whose bounding box
652 : // contain the point p.
653 406627 : for (auto c : index_range(children))
654 440397 : if (children[c]->bounds_point(p, relative_tol))
655 : {
656 15116 : const Elem * e =
657 92128 : children[c]->find_element(p,allowed_subdomains,
658 : relative_tol);
659 :
660 92128 : if (e != nullptr)
661 7556 : return e;
662 :
663 : // If we get here then a child that bounds the
664 : // point does not have any elements that contain
665 : // the point. So, we will search all our children.
666 : // However, we have already searched child c so there
667 : // is no use searching her again.
668 10 : searched_child[c] = true;
669 : }
670 :
671 :
672 : // If we get here then our child whose bounding box
673 : // was searched and did not find any elements containing
674 : // the point p. So, let's look at the other children
675 : // but exclude the one we have already searched.
676 50 : for (auto c : index_range(children))
677 40 : if (!searched_child[c])
678 : {
679 12 : const Elem * e =
680 30 : children[c]->find_element(p,allowed_subdomains,
681 : relative_tol);
682 :
683 30 : if (e != nullptr)
684 0 : return e;
685 : }
686 :
687 : // If we get here we have searched all our children.
688 : // Since this process was started at the root node then
689 : // we have searched all the elements in the tree without
690 : // success. So, we should return nullptr since at this point
691 : // _no_ elements in the tree claim to contain point p.
692 4 : return nullptr;
693 : }
694 :
695 :
696 :
697 : template <unsigned int N>
698 107916 : void TreeNode<N>::find_elements_in_children (const Point & p,
699 : std::set<const Elem *> & candidate_elements,
700 : const std::set<subdomain_id_type> * allowed_subdomains,
701 : Real relative_tol) const
702 : {
703 43334 : libmesh_assert (!this->active());
704 :
705 : // First only look in the children whose bounding box
706 : // contain the point p.
707 632768 : for (std::size_t c=0; c<children.size(); c++)
708 491396 : if (children[c]->bounds_point(p, relative_tol))
709 145180 : children[c]->find_elements(p, candidate_elements,
710 : allowed_subdomains, relative_tol);
711 107916 : }
712 :
713 :
714 :
715 : // ------------------------------------------------------------
716 : // Explicit Instantiations
717 : template class LIBMESH_EXPORT TreeNode<2>;
718 : template class LIBMESH_EXPORT TreeNode<4>;
719 : template class LIBMESH_EXPORT TreeNode<8>;
720 :
721 : } // namespace libMesh
|
