Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : #ifndef LIBMESH_ELEM_INTERNAL_H
19 : #define LIBMESH_ELEM_INTERNAL_H
20 :
21 : // libMesh includes
22 : #include "libmesh/elem.h"
23 :
24 : namespace libMesh
25 : {
26 :
27 : // Forward declarations
28 : class PointLocatorBase;
29 : class PeriodicBoundaries;
30 :
31 : /**
32 : * The ElemInternal namespace holds helper functions that are used
33 : * internally by the Elem class. These should not be called directly,
34 : * call the appropriate member functions on the Elem class instead.
35 : */
36 : namespace ElemInternal
37 : {
38 :
39 : template<class T>
40 : void
41 76113 : family_tree (T elem,
42 : std::vector<T> & family,
43 : bool reset = true)
44 : {
45 : // The "family tree" doesn't include subactive elements
46 6339 : libmesh_assert(!elem->subactive());
47 :
48 : // Clear the vector if the flag reset tells us to.
49 76113 : if (reset)
50 770 : family.clear();
51 :
52 : // Add this element to the family tree.
53 76113 : family.push_back(elem);
54 :
55 : // Recurse into the elements children, if it has them.
56 : // Do not clear the vector any more.
57 76113 : if (!elem->active())
58 75333 : for (auto & c : elem->child_ref_range())
59 66612 : if (!c.is_remote())
60 66612 : ElemInternal::family_tree (&c, family, false);
61 76113 : }
62 :
63 :
64 :
65 : template<class T>
66 : void
67 32308963 : total_family_tree(T elem,
68 : std::vector<T> & family,
69 : bool reset)
70 : {
71 : // Clear the vector if the flag reset tells us to.
72 32308963 : if (reset)
73 1111492 : family.clear();
74 :
75 : // Add this element to the family tree.
76 32308963 : family.push_back(elem);
77 :
78 : // Recurse into the elements children, if it has them.
79 : // Do not clear the vector any more.
80 32308963 : if (elem->has_children())
81 1847564 : for (auto & c : elem->child_ref_range())
82 1486120 : if (!c.is_remote())
83 1456333 : ElemInternal::total_family_tree (&c, family, false);
84 32308963 : }
85 :
86 :
87 :
88 : template<class T>
89 : void
90 11061623 : active_family_tree(T elem,
91 : std::vector<T> & active_family,
92 : bool reset = true)
93 : {
94 : // The "family tree" doesn't include subactive elements
95 292485 : libmesh_assert(!elem->subactive());
96 :
97 : // Clear the vector if the flag reset tells us to.
98 11061623 : if (reset)
99 45882 : active_family.clear();
100 :
101 : // Add this element to the family tree if it is active
102 11061623 : if (elem->active())
103 8945119 : active_family.push_back(elem);
104 :
105 : // Otherwise recurse into the element's children.
106 : // Do not clear the vector any more.
107 : else
108 12853160 : for (auto & c : elem->child_ref_range())
109 10736656 : if (!c.is_remote())
110 10027899 : ElemInternal::active_family_tree (&c, active_family, false);
111 11061623 : }
112 :
113 :
114 :
115 : template <class T>
116 : void
117 0 : family_tree_by_side (T elem,
118 : std::vector<T> & family,
119 : unsigned int s,
120 : bool reset)
121 : {
122 : // The "family tree" doesn't include subactive elements
123 0 : libmesh_assert(!elem->subactive());
124 :
125 : // Clear the vector if the flag reset tells us to.
126 0 : if (reset)
127 0 : family.clear();
128 :
129 0 : libmesh_assert_less (s, elem->n_sides());
130 :
131 : // Add this element to the family tree.
132 0 : family.push_back(elem);
133 :
134 : // Recurse into the elements children, if it has them.
135 : // Do not clear the vector any more.
136 0 : if (!elem->active())
137 : {
138 0 : const unsigned int nc = elem->n_children();
139 0 : for (unsigned int c = 0; c != nc; c++)
140 0 : if (!elem->child_ptr(c)->is_remote() && elem->is_child_on_side(c, s))
141 0 : ElemInternal::family_tree_by_side (elem->child_ptr(c), family, s, false);
142 : }
143 0 : }
144 :
145 :
146 :
147 : template <class T>
148 : void
149 284686783 : total_family_tree_by_side (T elem,
150 : std::vector<T> & family,
151 : unsigned int s,
152 : bool reset)
153 : {
154 : // Clear the vector if the flag reset tells us to.
155 284686783 : if (reset)
156 38682 : family.clear();
157 :
158 51062 : libmesh_assert_less (s, elem->n_sides());
159 :
160 : // Add this element to the family tree.
161 284686783 : family.push_back(elem);
162 :
163 : // Recurse into the elements children, if it has them.
164 : // Do not clear the vector any more.
165 284686783 : if (elem->has_children())
166 : {
167 52504007 : const unsigned int nc = elem->n_children();
168 293869715 : for (unsigned int c = 0; c != nc; c++)
169 241393062 : if (!elem->child_ptr(c)->is_remote() && elem->is_child_on_side(c, s))
170 64100639 : ElemInternal::total_family_tree_by_side (elem->child_ptr(c), family, s, false);
171 : }
172 284686783 : }
173 :
174 :
175 :
176 : template <class T>
177 : void
178 824200 : active_family_tree_by_side (T elem,
179 : std::vector<T> & family,
180 : unsigned int side,
181 : bool reset = true)
182 : {
183 : // The "family tree" doesn't include subactive or remote elements
184 66732 : libmesh_assert(!elem->subactive());
185 66732 : libmesh_assert(!elem->is_remote());
186 :
187 : // Clear the vector if the flag reset tells us to.
188 824200 : if (reset)
189 3045 : family.clear();
190 :
191 66732 : libmesh_assert_less (side, elem->n_sides());
192 :
193 : // Add an active element to the family tree.
194 824200 : if (elem->active())
195 745033 : family.push_back(elem);
196 :
197 : // Or recurse into an ancestor element's children.
198 : // Do not clear the vector any more.
199 : else
200 : {
201 79167 : const unsigned int nc = elem->n_children();
202 570577 : for (unsigned int c = 0; c != nc; c++)
203 530232 : if (!elem->child_ptr(c)->is_remote() && elem->is_child_on_side(c, side))
204 252488 : ElemInternal::active_family_tree_by_side(elem->child_ptr(c), family, side, false);
205 : }
206 824200 : }
207 :
208 :
209 : template <class T>
210 : void
211 0 : family_tree_by_neighbor (T elem,
212 : std::vector<T> & family,
213 : T neighbor_in,
214 : bool reset = true)
215 : {
216 : // The "family tree" doesn't include subactive elements
217 0 : libmesh_assert(!elem->subactive());
218 :
219 : // Clear the vector if the flag reset tells us to.
220 0 : if (reset)
221 0 : family.clear();
222 :
223 : // This only makes sense if we're already a neighbor
224 0 : libmesh_assert (elem->has_neighbor(neighbor_in));
225 :
226 : // Add this element to the family tree.
227 0 : family.push_back(elem);
228 :
229 : // Recurse into the elements children, if it's not active.
230 : // Do not clear the vector any more.
231 0 : if (!elem->active())
232 0 : for (auto & c : elem->child_ref_range())
233 0 : if (!c.is_remote() && c.has_neighbor(neighbor_in))
234 0 : ElemInternal::family_tree_by_neighbor (&c, family, neighbor_in, false);
235 0 : }
236 :
237 :
238 : template <class T>
239 : void
240 70152520 : total_family_tree_by_neighbor (T elem,
241 : std::vector<T> & family,
242 : T neighbor_in,
243 : bool reset = true)
244 : {
245 : // Clear the vector if the flag reset tells us to.
246 70152520 : if (reset)
247 3917 : family.clear();
248 :
249 : // This only makes sense if we're already a neighbor
250 3917 : libmesh_assert (elem->has_neighbor(neighbor_in));
251 :
252 : // Add this element to the family tree.
253 70152520 : family.push_back(elem);
254 :
255 : // Recurse into the elements children, if it has any.
256 : // Do not clear the vector any more.
257 70152520 : if (elem->has_children())
258 61178660 : for (auto & c : elem->child_ref_range())
259 49264466 : if (!c.is_remote() && c.has_neighbor(neighbor_in))
260 7787 : ElemInternal::total_family_tree_by_neighbor (&c, family, neighbor_in, false);
261 70152520 : }
262 :
263 :
264 : template <class T>
265 : void
266 0 : family_tree_by_subneighbor (T elem,
267 : std::vector<T> & family,
268 : T neighbor_in,
269 : T subneighbor,
270 : bool reset = true)
271 : {
272 : // The "family tree" doesn't include subactive elements
273 0 : libmesh_assert(!elem->subactive());
274 :
275 : // Clear the vector if the flag reset tells us to.
276 0 : if (reset)
277 0 : family.clear();
278 :
279 : // To simplify this function we need an existing neighbor
280 0 : libmesh_assert (neighbor_in);
281 0 : libmesh_assert (!neighbor_in->is_remote());
282 0 : libmesh_assert (elem->has_neighbor(neighbor_in));
283 :
284 : // This only makes sense if subneighbor descends from neighbor
285 0 : libmesh_assert (subneighbor);
286 0 : libmesh_assert (!subneighbor->is_remote());
287 0 : libmesh_assert (neighbor_in->is_ancestor_of(subneighbor));
288 :
289 : // Add this element to the family tree if applicable.
290 0 : if (neighbor_in == subneighbor)
291 0 : family.push_back(elem);
292 :
293 : // Recurse into the elements children, if it's not active.
294 : // Do not clear the vector any more.
295 0 : if (!elem->active())
296 0 : for (auto & c : elem->child_ref_range())
297 0 : if (!c.is_remote())
298 0 : for (auto child_neigh : c.neighbor_ptr_range())
299 0 : if (child_neigh &&
300 0 : (child_neigh == neighbor_in || (child_neigh->parent() == neighbor_in &&
301 : child_neigh->is_ancestor_of(subneighbor))))
302 0 : ElemInternal::family_tree_by_subneighbor (&c, family, child_neigh, subneighbor, false);
303 0 : }
304 :
305 :
306 :
307 : template <class T>
308 : void
309 2874 : total_family_tree_by_subneighbor(T elem,
310 : std::vector<T> & family,
311 : T neighbor_in,
312 : T subneighbor,
313 : bool reset = true)
314 : {
315 : // Clear the vector if the flag reset tells us to.
316 2874 : if (reset)
317 0 : family.clear();
318 :
319 : // To simplify this function we need an existing neighbor
320 0 : libmesh_assert (neighbor_in);
321 0 : libmesh_assert (!neighbor_in->is_remote());
322 0 : libmesh_assert (elem->has_neighbor(neighbor_in));
323 :
324 : // This only makes sense if subneighbor descends from neighbor
325 0 : libmesh_assert (subneighbor);
326 0 : libmesh_assert (!subneighbor->is_remote());
327 0 : libmesh_assert (neighbor_in->is_ancestor_of(subneighbor));
328 :
329 : // Add this element to the family tree if applicable.
330 2874 : if (neighbor_in == subneighbor)
331 1114 : family.push_back(elem);
332 :
333 : // Recurse into the elements children, if it has any.
334 : // Do not clear the vector any more.
335 2874 : if (elem->has_children())
336 12740 : for (auto & c : elem->child_ref_range())
337 10980 : if (!c.is_remote())
338 51274 : for (auto child_neigh : c.neighbor_ptr_range())
339 43618 : if (child_neigh &&
340 41585 : (child_neigh == neighbor_in ||
341 3102 : (child_neigh->parent() == neighbor_in && child_neigh->is_ancestor_of(subneighbor))))
342 1114 : c.total_family_tree_by_subneighbor(family, child_neigh, subneighbor, false);
343 2874 : }
344 :
345 :
346 :
347 : template <class T>
348 : void
349 51481142 : active_family_tree_by_neighbor(T elem,
350 : std::vector<T> & family,
351 : T neighbor_in,
352 : bool reset = true)
353 : {
354 : // The "family tree" doesn't include subactive elements or
355 : // remote_elements
356 2987596 : libmesh_assert(!elem->subactive());
357 2987596 : libmesh_assert(!elem->is_remote());
358 :
359 : // Clear the vector if the flag reset tells us to.
360 51481142 : if (reset)
361 2696583 : family.clear();
362 :
363 : // This only makes sense if we're already a neighbor
364 : #ifndef NDEBUG
365 2987596 : if (elem->level() >= neighbor_in->level())
366 2841523 : libmesh_assert (elem->has_neighbor(neighbor_in));
367 : #endif
368 :
369 : // Add an active element to the family tree.
370 51481142 : if (elem->active())
371 48254731 : family.push_back(elem);
372 :
373 : // Or recurse into an ancestor element's children.
374 : // Do not clear the vector any more.
375 142021 : else if (!elem->active())
376 21980103 : for (auto & c : elem->child_ref_range())
377 19335402 : if (!c.is_remote() && c.has_neighbor(neighbor_in))
378 9231972 : ElemInternal::active_family_tree_by_neighbor (&c, family, neighbor_in, false);
379 51481142 : }
380 :
381 :
382 :
383 : template <class T>
384 : void
385 26334 : active_family_tree_by_topological_neighbor (T elem,
386 : std::vector<T> & family,
387 : T neighbor_in,
388 : const MeshBase & mesh,
389 : const PointLocatorBase & point_locator,
390 : const PeriodicBoundaries * pb,
391 : bool reset = true)
392 : {
393 : // The "family tree" doesn't include subactive elements or
394 : // remote_elements
395 10514 : libmesh_assert(!elem->subactive());
396 10514 : libmesh_assert(!elem->is_remote());
397 :
398 : // Clear the vector if the flag reset tells us to.
399 26334 : if (reset)
400 5290 : family.clear();
401 :
402 : // This only makes sense if we're already a topological neighbor
403 : #ifndef NDEBUG
404 10514 : if (elem->level() >= neighbor_in->level())
405 8718 : libmesh_assert (elem->has_topological_neighbor(neighbor_in,
406 : mesh,
407 : point_locator,
408 : pb));
409 : #endif
410 :
411 : // Add an active element to the family tree.
412 26334 : if (elem->active())
413 23512 : family.push_back(elem);
414 :
415 : // Or recurse into an ancestor element's children.
416 : // Do not clear the vector any more.
417 1306 : else if (!elem->active())
418 14110 : for (auto & c : elem->child_ref_range())
419 12128 : if (!c.is_remote() &&
420 840 : c.has_topological_neighbor(neighbor_in, mesh, point_locator, pb))
421 10448 : c.active_family_tree_by_topological_neighbor
422 840 : (family, neighbor_in, mesh, point_locator, pb, false);
423 26334 : }
424 :
425 :
426 :
427 : template <class T>
428 : void
429 9823332 : find_point_neighbors(T this_elem,
430 : std::set<T> & neighbor_set,
431 : T start_elem)
432 : {
433 16767 : libmesh_assert(start_elem);
434 16767 : libmesh_assert(start_elem->active());
435 16767 : libmesh_assert(start_elem->contains_vertex_of(this_elem, true) ||
436 : this_elem->contains_vertex_of(start_elem, true));
437 :
438 16767 : neighbor_set.clear();
439 9806565 : neighbor_set.insert(start_elem);
440 :
441 33534 : std::set<T> untested_set, next_untested_set;
442 9806565 : untested_set.insert(start_elem);
443 :
444 60506569 : while (!untested_set.empty())
445 : {
446 : // Loop over all the elements in the patch that haven't already
447 : // been tested
448 333344684 : for (const auto & elem : untested_set)
449 1466720573 : for (auto current_neighbor : elem->neighbor_ptr_range())
450 : {
451 2316570367 : if (current_neighbor &&
452 1132903825 : !current_neighbor->is_remote()) // we have a real neighbor on this side
453 : {
454 3392790 : if (current_neighbor->active()) // ... if it is active
455 : {
456 1095590616 : if (this_elem->contains_vertex_of(current_neighbor, true) // ... and touches us
457 1097286783 : || current_neighbor->contains_vertex_of(this_elem, true))
458 : {
459 : // Make sure we'll test it
460 1293470 : if (!neighbor_set.count(current_neighbor))
461 271609568 : next_untested_set.insert (current_neighbor);
462 :
463 : // And add it
464 875139960 : neighbor_set.insert (current_neighbor);
465 : }
466 : }
467 : #ifdef LIBMESH_ENABLE_AMR
468 : else // ... the neighbor is *not* active,
469 : { // ... so add *all* neighboring
470 : // active children
471 456 : std::vector<T> active_neighbor_children;
472 :
473 456 : current_neighbor->active_family_tree_by_neighbor
474 2116517 : (active_neighbor_children, elem);
475 :
476 9021465 : for (const auto & current_child : active_neighbor_children)
477 : {
478 10959184 : if (this_elem->contains_vertex_of(current_child, true) ||
479 4054692 : current_child->contains_vertex_of(this_elem, true))
480 : {
481 : // Make sure we'll test it
482 328 : if (!neighbor_set.count(current_child))
483 852730 : next_untested_set.insert (current_child);
484 :
485 2848704 : neighbor_set.insert (current_child);
486 : }
487 : }
488 : }
489 : #endif // #ifdef LIBMESH_ENABLE_AMR
490 : }
491 : }
492 81353 : untested_set.swap(next_untested_set);
493 81353 : next_untested_set.clear();
494 : }
495 9823332 : }
496 :
497 :
498 :
499 : template <class T>
500 : void
501 11629 : find_interior_neighbors(T this_elem,
502 : std::set<T> & neighbor_set)
503 : {
504 1268 : neighbor_set.clear();
505 :
506 23258 : if ((this_elem->dim() >= LIBMESH_DIM) ||
507 11629 : !this_elem->interior_parent())
508 0 : return;
509 :
510 11629 : T ip = this_elem->interior_parent();
511 1268 : libmesh_assert (ip->contains_vertex_of(this_elem, true) ||
512 : this_elem->contains_vertex_of(ip, true));
513 :
514 1268 : libmesh_assert (!ip->subactive());
515 :
516 : #ifdef LIBMESH_ENABLE_AMR
517 1268 : while (!ip->active()) // only possible with AMR, be careful because
518 : { // ip->child_ptr(c) is only good with AMR.
519 0 : for (auto & child : ip->child_ref_range())
520 : {
521 0 : if (child.contains_vertex_of(this_elem, true) ||
522 0 : this_elem->contains_vertex_of(&child, true))
523 : {
524 0 : ip = &child;
525 0 : break;
526 : }
527 : }
528 : }
529 : #endif
530 :
531 11629 : ElemInternal::find_point_neighbors(this_elem, neighbor_set, ip);
532 :
533 : // Now we have all point neighbors from the interior manifold, but
534 : // we need to weed out any neighbors that *only* intersect us at one
535 : // point (or at one edge, if we're a 1-D element in 3D).
536 : //
537 : // The refinement hierarchy helps us here: if the interior element
538 : // has a lower or equal refinement level then we can discard it iff
539 : // it doesn't contain all our vertices. If it has a higher
540 : // refinement level then we can discard it iff we don't contain at
541 : // least dim()+1 of its vertices
542 1268 : auto it = neighbor_set.begin();
543 1268 : const auto end = neighbor_set.end();
544 :
545 68201 : while (it != end)
546 : {
547 4524 : auto current = it++;
548 56572 : T current_elem = *current;
549 :
550 : // This won't invalidate iterator it, because it is already
551 : // pointing to the next element
552 56572 : if (current_elem->level() > this_elem->level())
553 : {
554 28 : unsigned int vertices_contained = 0;
555 1940 : for (auto & n : current_elem->node_ref_range())
556 1746 : if (this_elem->contains_point(n))
557 194 : vertices_contained++;
558 :
559 194 : if (vertices_contained <= this_elem->dim())
560 : {
561 166 : neighbor_set.erase(current);
562 28 : continue;
563 : }
564 : }
565 : else
566 : {
567 123024 : for (auto & n : this_elem->node_ref_range())
568 : {
569 106790 : if (!current_elem->contains_point(n))
570 : {
571 37092 : neighbor_set.erase(current);
572 3052 : break;
573 : }
574 : }
575 : }
576 : }
577 : }
578 :
579 : } // namespace ElemInternal
580 : } // namespace libMesh
581 :
582 : #endif
|
