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 : #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 29196678 : 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 29196678 : if (reset)
73 1076666 : family.clear();
74 :
75 : // Add this element to the family tree.
76 29196678 : family.push_back(elem);
77 :
78 : // Recurse into the elements children, if it has them.
79 : // Do not clear the vector any more.
80 29196678 : if (elem->has_children())
81 1597904 : for (auto & c : elem->child_ref_range())
82 1286392 : if (!c.is_remote())
83 1256605 : ElemInternal::total_family_tree (&c, family, false);
84 29196678 : }
85 :
86 :
87 :
88 : template<class T>
89 : void
90 10358061 : 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 247197 : libmesh_assert(!elem->subactive());
96 :
97 : // Clear the vector if the flag reset tells us to.
98 10358061 : if (reset)
99 38110 : active_family.clear();
100 :
101 : // Add this element to the family tree if it is active
102 10358061 : if (elem->active())
103 8378699 : 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 12023066 : for (auto & c : elem->child_ref_range())
109 10043704 : if (!c.is_remote())
110 9431011 : ElemInternal::active_family_tree (&c, active_family, false);
111 10358061 : }
112 :
113 :
114 :
115 : template <class T>
116 : void
117 41710533 : 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 638 : libmesh_assert(!elem->subactive());
124 :
125 : // Clear the vector if the flag reset tells us to.
126 41710533 : if (reset)
127 638 : family.clear();
128 :
129 638 : libmesh_assert_less (s, elem->n_sides());
130 :
131 : // Add this element to the family tree.
132 41710533 : family.push_back(elem);
133 :
134 : // Recurse into the elements children, if it has them.
135 : // Do not clear the vector any more.
136 41710533 : if (!elem->active())
137 : {
138 138733 : const unsigned int nc = elem->n_children();
139 883057 : for (unsigned int c = 0; c != nc; c++)
140 744324 : if (!elem->child_ptr(c)->is_remote() && elem->is_child_on_side(c, s))
141 12449 : ElemInternal::family_tree_by_side (elem->child_ptr(c), family, s, false);
142 : }
143 41710533 : }
144 :
145 :
146 :
147 : template <class T>
148 : void
149 829214 : active_family_tree_by_side (T elem,
150 : std::vector<T> & family,
151 : unsigned int side,
152 : bool reset = true)
153 : {
154 : // The "family tree" doesn't include subactive or remote elements
155 69594 : libmesh_assert(!elem->subactive());
156 69594 : libmesh_assert(!elem->is_remote());
157 :
158 : // Clear the vector if the flag reset tells us to.
159 829214 : if (reset)
160 3037 : family.clear();
161 :
162 69594 : libmesh_assert_less (side, elem->n_sides());
163 :
164 : // Add an active element to the family tree.
165 829214 : if (elem->active())
166 750047 : family.push_back(elem);
167 :
168 : // Or recurse into an ancestor element's children.
169 : // Do not clear the vector any more.
170 : else
171 : {
172 79167 : const unsigned int nc = elem->n_children();
173 570577 : for (unsigned int c = 0; c != nc; c++)
174 530232 : if (!elem->child_ptr(c)->is_remote() && elem->is_child_on_side(c, side))
175 252488 : ElemInternal::active_family_tree_by_side(elem->child_ptr(c), family, side, false);
176 : }
177 829214 : }
178 :
179 :
180 : template <class T>
181 : void
182 0 : family_tree_by_neighbor (T elem,
183 : std::vector<T> & family,
184 : T neighbor_in,
185 : bool reset = true)
186 : {
187 : // The "family tree" doesn't include subactive elements
188 0 : libmesh_assert(!elem->subactive());
189 :
190 : // Clear the vector if the flag reset tells us to.
191 0 : if (reset)
192 0 : family.clear();
193 :
194 : // This only makes sense if we're already a neighbor
195 0 : libmesh_assert (elem->has_neighbor(neighbor_in));
196 :
197 : // Add this element to the family tree.
198 0 : family.push_back(elem);
199 :
200 : // Recurse into the elements children, if it's not active.
201 : // Do not clear the vector any more.
202 0 : if (!elem->active())
203 0 : for (auto & c : elem->child_ref_range())
204 0 : if (!c.is_remote() && c.has_neighbor(neighbor_in))
205 0 : ElemInternal::family_tree_by_neighbor (&c, family, neighbor_in, false);
206 0 : }
207 :
208 :
209 : template <class T>
210 : void
211 68505719 : total_family_tree_by_neighbor (T elem,
212 : std::vector<T> & family,
213 : T neighbor_in,
214 : bool reset = true)
215 : {
216 : // Clear the vector if the flag reset tells us to.
217 68505719 : if (reset)
218 3619 : family.clear();
219 :
220 : // This only makes sense if we're already a neighbor
221 3619 : libmesh_assert (elem->has_neighbor(neighbor_in));
222 :
223 : // Add this element to the family tree.
224 68505719 : family.push_back(elem);
225 :
226 : // Recurse into the elements children, if it has any.
227 : // Do not clear the vector any more.
228 68505719 : if (elem->has_children())
229 61082221 : for (auto & c : elem->child_ref_range())
230 49183486 : if (!c.is_remote() && c.has_neighbor(neighbor_in))
231 7783 : ElemInternal::total_family_tree_by_neighbor (&c, family, neighbor_in, false);
232 68505719 : }
233 :
234 :
235 : template <class T>
236 : void
237 0 : family_tree_by_subneighbor (T elem,
238 : std::vector<T> & family,
239 : T neighbor_in,
240 : T subneighbor,
241 : bool reset = true)
242 : {
243 : // The "family tree" doesn't include subactive elements
244 0 : libmesh_assert(!elem->subactive());
245 :
246 : // Clear the vector if the flag reset tells us to.
247 0 : if (reset)
248 0 : family.clear();
249 :
250 : // To simplify this function we need an existing neighbor
251 0 : libmesh_assert (neighbor_in);
252 0 : libmesh_assert (!neighbor_in->is_remote());
253 0 : libmesh_assert (elem->has_neighbor(neighbor_in));
254 :
255 : // This only makes sense if subneighbor descends from neighbor
256 0 : libmesh_assert (subneighbor);
257 0 : libmesh_assert (!subneighbor->is_remote());
258 0 : libmesh_assert (neighbor_in->is_ancestor_of(subneighbor));
259 :
260 : // Add this element to the family tree if applicable.
261 0 : if (neighbor_in == subneighbor)
262 0 : family.push_back(elem);
263 :
264 : // Recurse into the elements children, if it's not active.
265 : // Do not clear the vector any more.
266 0 : if (!elem->active())
267 0 : for (auto & c : elem->child_ref_range())
268 0 : if (!c.is_remote())
269 0 : for (auto child_neigh : c.neighbor_ptr_range())
270 0 : if (child_neigh &&
271 0 : (child_neigh == neighbor_in || (child_neigh->parent() == neighbor_in &&
272 : child_neigh->is_ancestor_of(subneighbor))))
273 0 : ElemInternal::family_tree_by_subneighbor (&c, family, child_neigh, subneighbor, false);
274 0 : }
275 :
276 :
277 :
278 : template <class T>
279 : void
280 2874 : total_family_tree_by_subneighbor(T elem,
281 : std::vector<T> & family,
282 : T neighbor_in,
283 : T subneighbor,
284 : bool reset = true)
285 : {
286 : // Clear the vector if the flag reset tells us to.
287 2874 : if (reset)
288 0 : family.clear();
289 :
290 : // To simplify this function we need an existing neighbor
291 0 : libmesh_assert (neighbor_in);
292 0 : libmesh_assert (!neighbor_in->is_remote());
293 0 : libmesh_assert (elem->has_neighbor(neighbor_in));
294 :
295 : // This only makes sense if subneighbor descends from neighbor
296 0 : libmesh_assert (subneighbor);
297 0 : libmesh_assert (!subneighbor->is_remote());
298 0 : libmesh_assert (neighbor_in->is_ancestor_of(subneighbor));
299 :
300 : // Add this element to the family tree if applicable.
301 2874 : if (neighbor_in == subneighbor)
302 1114 : family.push_back(elem);
303 :
304 : // Recurse into the elements children, if it has any.
305 : // Do not clear the vector any more.
306 2874 : if (elem->has_children())
307 12740 : for (auto & c : elem->child_ref_range())
308 10980 : if (!c.is_remote())
309 51274 : for (auto child_neigh : c.neighbor_ptr_range())
310 43618 : if (child_neigh &&
311 41585 : (child_neigh == neighbor_in ||
312 3102 : (child_neigh->parent() == neighbor_in && child_neigh->is_ancestor_of(subneighbor))))
313 1114 : c.total_family_tree_by_subneighbor(family, child_neigh, subneighbor, false);
314 2874 : }
315 :
316 :
317 :
318 : template <class T>
319 : void
320 48042327 : active_family_tree_by_neighbor(T elem,
321 : std::vector<T> & family,
322 : T neighbor_in,
323 : bool reset = true)
324 : {
325 : // The "family tree" doesn't include subactive elements or
326 : // remote_elements
327 2684359 : libmesh_assert(!elem->subactive());
328 2684359 : libmesh_assert(!elem->is_remote());
329 :
330 : // Clear the vector if the flag reset tells us to.
331 48042327 : if (reset)
332 2438638 : family.clear();
333 :
334 : // This only makes sense if we're already a neighbor
335 : #ifndef NDEBUG
336 2684359 : if (elem->level() >= neighbor_in->level())
337 2583834 : libmesh_assert (elem->has_neighbor(neighbor_in));
338 : #endif
339 :
340 : // Add an active element to the family tree.
341 48042327 : if (elem->active())
342 44991532 : family.push_back(elem);
343 :
344 : // Or recurse into an ancestor element's children.
345 : // Do not clear the vector any more.
346 119374 : else if (!elem->active())
347 21130249 : for (auto & c : elem->child_ref_range())
348 18570812 : if (!c.is_remote() && c.has_neighbor(neighbor_in))
349 8881273 : ElemInternal::active_family_tree_by_neighbor (&c, family, neighbor_in, false);
350 48042327 : }
351 :
352 :
353 :
354 : template <class T>
355 : void
356 26334 : active_family_tree_by_topological_neighbor (T elem,
357 : std::vector<T> & family,
358 : T neighbor_in,
359 : const MeshBase & mesh,
360 : const PointLocatorBase & point_locator,
361 : const PeriodicBoundaries * pb,
362 : bool reset = true)
363 : {
364 : // The "family tree" doesn't include subactive elements or
365 : // remote_elements
366 10514 : libmesh_assert(!elem->subactive());
367 10514 : libmesh_assert(!elem->is_remote());
368 :
369 : // Clear the vector if the flag reset tells us to.
370 26334 : if (reset)
371 5290 : family.clear();
372 :
373 : // This only makes sense if we're already a topological neighbor
374 : #ifndef NDEBUG
375 10514 : if (elem->level() >= neighbor_in->level())
376 8718 : libmesh_assert (elem->has_topological_neighbor(neighbor_in,
377 : mesh,
378 : point_locator,
379 : pb));
380 : #endif
381 :
382 : // Add an active element to the family tree.
383 26334 : if (elem->active())
384 23512 : family.push_back(elem);
385 :
386 : // Or recurse into an ancestor element's children.
387 : // Do not clear the vector any more.
388 1306 : else if (!elem->active())
389 14110 : for (auto & c : elem->child_ref_range())
390 12128 : if (!c.is_remote() &&
391 840 : c.has_topological_neighbor(neighbor_in, mesh, point_locator, pb))
392 10448 : c.active_family_tree_by_topological_neighbor
393 840 : (family, neighbor_in, mesh, point_locator, pb, false);
394 26334 : }
395 :
396 :
397 :
398 : template <class T>
399 : void
400 9419078 : find_point_neighbors(T this_elem,
401 : std::set<T> & neighbor_set,
402 : T start_elem)
403 : {
404 15149 : libmesh_assert(start_elem);
405 15149 : libmesh_assert(start_elem->active());
406 15149 : libmesh_assert(start_elem->contains_vertex_of(this_elem, true) ||
407 : this_elem->contains_vertex_of(start_elem, true));
408 :
409 15149 : neighbor_set.clear();
410 9403929 : neighbor_set.insert(start_elem);
411 :
412 30298 : std::set<T> untested_set, next_untested_set;
413 9403929 : untested_set.insert(start_elem);
414 :
415 57379359 : while (!untested_set.empty())
416 : {
417 : // Loop over all the elements in the patch that haven't already
418 : // been tested
419 316224370 : for (const auto & elem : untested_set)
420 1386288077 : for (auto current_neighbor : elem->neighbor_ptr_range())
421 : {
422 2188055034 : if (current_neighbor &&
423 1070303882 : !current_neighbor->is_remote()) // we have a real neighbor on this side
424 : {
425 3071746 : if (current_neighbor->active()) // ... if it is active
426 : {
427 1036782959 : if (this_elem->contains_vertex_of(current_neighbor, true) // ... and touches us
428 1038318604 : || current_neighbor->contains_vertex_of(this_elem, true))
429 : {
430 : // Make sure we'll test it
431 1182546 : if (!neighbor_set.count(current_neighbor))
432 257538395 : next_untested_set.insert (current_neighbor);
433 :
434 : // And add it
435 828471850 : neighbor_set.insert (current_neighbor);
436 : }
437 : }
438 : #ifdef LIBMESH_ENABLE_AMR
439 : else // ... the neighbor is *not* active,
440 : { // ... so add *all* neighboring
441 : // active children
442 456 : std::vector<T> active_neighbor_children;
443 :
444 456 : current_neighbor->active_family_tree_by_neighbor
445 2150099 : (active_neighbor_children, elem);
446 :
447 9122964 : for (const auto & current_child : active_neighbor_children)
448 : {
449 11046574 : if (this_elem->contains_vertex_of(current_child, true) ||
450 4074165 : current_child->contains_vertex_of(this_elem, true))
451 : {
452 : // Make sure we'll test it
453 328 : if (!neighbor_set.count(current_child))
454 867801 : next_untested_set.insert (current_child);
455 :
456 2897148 : neighbor_set.insert (current_child);
457 : }
458 : }
459 : }
460 : #endif // #ifdef LIBMESH_ENABLE_AMR
461 : }
462 : }
463 75415 : untested_set.swap(next_untested_set);
464 75415 : next_untested_set.clear();
465 : }
466 9419078 : }
467 :
468 :
469 :
470 : template <class T>
471 : void
472 11629 : find_interior_neighbors(T this_elem,
473 : std::set<T> & neighbor_set)
474 : {
475 1268 : neighbor_set.clear();
476 :
477 23258 : if ((this_elem->dim() >= LIBMESH_DIM) ||
478 11629 : !this_elem->interior_parent())
479 0 : return;
480 :
481 11629 : T ip = this_elem->interior_parent();
482 1268 : libmesh_assert (ip->contains_vertex_of(this_elem, true) ||
483 : this_elem->contains_vertex_of(ip, true));
484 :
485 1268 : libmesh_assert (!ip->subactive());
486 :
487 : #ifdef LIBMESH_ENABLE_AMR
488 1268 : while (!ip->active()) // only possible with AMR, be careful because
489 : { // ip->child_ptr(c) is only good with AMR.
490 0 : for (auto & child : ip->child_ref_range())
491 : {
492 0 : if (child.contains_vertex_of(this_elem, true) ||
493 0 : this_elem->contains_vertex_of(&child, true))
494 : {
495 0 : ip = &child;
496 0 : break;
497 : }
498 : }
499 : }
500 : #endif
501 :
502 11629 : ElemInternal::find_point_neighbors(this_elem, neighbor_set, ip);
503 :
504 : // Now we have all point neighbors from the interior manifold, but
505 : // we need to weed out any neighbors that *only* intersect us at one
506 : // point (or at one edge, if we're a 1-D element in 3D).
507 : //
508 : // The refinement hierarchy helps us here: if the interior element
509 : // has a lower or equal refinement level then we can discard it iff
510 : // it doesn't contain all our vertices. If it has a higher
511 : // refinement level then we can discard it iff we don't contain at
512 : // least dim()+1 of its vertices
513 1268 : auto it = neighbor_set.begin();
514 1268 : const auto end = neighbor_set.end();
515 :
516 68201 : while (it != end)
517 : {
518 4524 : auto current = it++;
519 56572 : T current_elem = *current;
520 :
521 : // This won't invalidate iterator it, because it is already
522 : // pointing to the next element
523 56572 : if (current_elem->level() > this_elem->level())
524 : {
525 28 : unsigned int vertices_contained = 0;
526 1940 : for (auto & n : current_elem->node_ref_range())
527 1746 : if (this_elem->contains_point(n))
528 194 : vertices_contained++;
529 :
530 194 : if (vertices_contained <= this_elem->dim())
531 : {
532 166 : neighbor_set.erase(current);
533 28 : continue;
534 : }
535 : }
536 : else
537 : {
538 123024 : for (auto & n : this_elem->node_ref_range())
539 : {
540 106790 : if (!current_elem->contains_point(n))
541 : {
542 37092 : neighbor_set.erase(current);
543 3052 : break;
544 : }
545 : }
546 : }
547 : }
548 : }
549 :
550 : } // namespace ElemInternal
551 : } // namespace libMesh
552 :
553 : #endif
|
