1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 family_tree (T elem,
42 std::vector<T> & family,
43 bool reset = true)
44 {
45 // The "family tree" doesn't include subactive elements
46 libmesh_assert(!elem->subactive());
47
48 // Clear the vector if the flag reset tells us to.
49 if (reset)
50 family.clear();
51
52 // Add this element to the family tree.
53 family.push_back(elem);
54
55 // Recurse into the elements children, if it has them.
56 // Do not clear the vector any more.
57 if (!elem->active())
58 for (auto & c : elem->child_ref_range())
59 if (!c.is_remote())
60 ElemInternal::family_tree (&c, family, false);
61 }
62
63
64
65 template<class T>
66 void
total_family_tree(T elem,std::vector<T> & family,bool reset)67 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 if (reset)
73 family.clear();
74
75 // Add this element to the family tree.
76 family.push_back(elem);
77
78 // Recurse into the elements children, if it has them.
79 // Do not clear the vector any more.
80 if (elem->has_children())
81 for (auto & c : elem->child_ref_range())
82 if (!c.is_remote())
83 ElemInternal::total_family_tree (&c, family, false);
84 }
85
86
87
88 template<class T>
89 void
90 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 libmesh_assert(!elem->subactive());
96
97 // Clear the vector if the flag reset tells us to.
98 if (reset)
99 active_family.clear();
100
101 // Add this element to the family tree if it is active
102 if (elem->active())
103 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 for (auto & c : elem->child_ref_range())
109 if (!c.is_remote())
110 ElemInternal::active_family_tree (&c, active_family, false);
111 }
112
113
114
115 template <class T>
116 void
family_tree_by_side(T elem,std::vector<T> & family,unsigned int s,bool reset)117 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 libmesh_assert(!elem->subactive());
124
125 // Clear the vector if the flag reset tells us to.
126 if (reset)
127 family.clear();
128
129 libmesh_assert_less (s, elem->n_sides());
130
131 // Add this element to the family tree.
132 family.push_back(elem);
133
134 // Recurse into the elements children, if it has them.
135 // Do not clear the vector any more.
136 if (!elem->active())
137 {
138 const unsigned int nc = elem->n_children();
139 for (unsigned int c = 0; c != nc; c++)
140 if (!elem->child_ptr(c)->is_remote() && elem->is_child_on_side(c, s))
141 ElemInternal::family_tree_by_side (elem->child_ptr(c), family, s, false);
142 }
143 }
144
145
146
147 template <class T>
148 void
149 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 libmesh_assert(!elem->subactive());
156 libmesh_assert(!elem->is_remote());
157
158 // Clear the vector if the flag reset tells us to.
159 if (reset)
160 family.clear();
161
162 libmesh_assert_less (side, elem->n_sides());
163
164 // Add an active element to the family tree.
165 if (elem->active())
166 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 const unsigned int nc = elem->n_children();
173 for (unsigned int c = 0; c != nc; c++)
174 if (!elem->child_ptr(c)->is_remote() && elem->is_child_on_side(c, side))
175 ElemInternal::active_family_tree_by_side(elem->child_ptr(c), family, side, false);
176 }
177 }
178
179
180 template <class T>
181 void
182 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 libmesh_assert(!elem->subactive());
189
190 // Clear the vector if the flag reset tells us to.
191 if (reset)
192 family.clear();
193
194 // This only makes sense if we're already a neighbor
195 libmesh_assert (elem->has_neighbor(neighbor_in));
196
197 // Add this element to the family tree.
198 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 if (!elem->active())
203 for (auto & c : elem->child_ref_range())
204 if (!c.is_remote() && c.has_neighbor(neighbor_in))
205 ElemInternal::family_tree_by_neighbor (&c, family, neighbor_in, false);
206 }
207
208
209 template <class T>
210 void
211 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 if (reset)
218 family.clear();
219
220 // This only makes sense if we're already a neighbor
221 libmesh_assert (elem->has_neighbor(neighbor_in));
222
223 // Add this element to the family tree.
224 family.push_back(elem);
225
226 // Recurse into the elements children, if it has any.
227 // Do not clear the vector any more.
228 if (elem->has_children())
229 for (auto & c : elem->child_ref_range())
230 if (!c.is_remote() && c.has_neighbor(neighbor_in))
231 ElemInternal::total_family_tree_by_neighbor (&c, family, neighbor_in, false);
232 }
233
234
235 template <class T>
236 void
237 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 libmesh_assert(!elem->subactive());
245
246 // Clear the vector if the flag reset tells us to.
247 if (reset)
248 family.clear();
249
250 // To simplify this function we need an existing neighbor
251 libmesh_assert (neighbor_in);
252 libmesh_assert (!neighbor_in->is_remote());
253 libmesh_assert (elem->has_neighbor(neighbor_in));
254
255 // This only makes sense if subneighbor descends from neighbor
256 libmesh_assert (subneighbor);
257 libmesh_assert (!subneighbor->is_remote());
258 libmesh_assert (neighbor_in->is_ancestor_of(subneighbor));
259
260 // Add this element to the family tree if applicable.
261 if (neighbor_in == subneighbor)
262 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 if (!elem->active())
267 for (auto & c : elem->child_ref_range())
268 if (!c.is_remote())
269 for (auto child_neigh : c.neighbor_ptr_range())
270 if (child_neigh &&
271 (child_neigh == neighbor_in || (child_neigh->parent() == neighbor_in &&
272 child_neigh->is_ancestor_of(subneighbor))))
273 ElemInternal::family_tree_by_subneighbor (&c, family, child_neigh, subneighbor, false);
274 }
275
276
277
278 template <class T>
279 void
280 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 if (reset)
288 family.clear();
289
290 // To simplify this function we need an existing neighbor
291 libmesh_assert (neighbor_in);
292 libmesh_assert (!neighbor_in->is_remote());
293 libmesh_assert (elem->has_neighbor(neighbor_in));
294
295 // This only makes sense if subneighbor descends from neighbor
296 libmesh_assert (subneighbor);
297 libmesh_assert (!subneighbor->is_remote());
298 libmesh_assert (neighbor_in->is_ancestor_of(subneighbor));
299
300 // Add this element to the family tree if applicable.
301 if (neighbor_in == subneighbor)
302 family.push_back(elem);
303
304 // Recurse into the elements children, if it has any.
305 // Do not clear the vector any more.
306 if (elem->has_children())
307 for (auto & c : elem->child_ref_range())
308 if (!c.is_remote())
309 for (auto child_neigh : c.neighbor_ptr_range())
310 if (child_neigh &&
311 (child_neigh == neighbor_in ||
312 (child_neigh->parent() == neighbor_in && child_neigh->is_ancestor_of(subneighbor))))
313 c.total_family_tree_by_subneighbor(family, child_neigh, subneighbor, false);
314 }
315
316
317
318 template <class T>
319 void
320 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 libmesh_assert(!elem->subactive());
328 libmesh_assert(!elem->is_remote());
329
330 // Clear the vector if the flag reset tells us to.
331 if (reset)
332 family.clear();
333
334 // This only makes sense if we're already a neighbor
335 #ifndef NDEBUG
336 if (elem->level() >= neighbor_in->level())
337 libmesh_assert (elem->has_neighbor(neighbor_in));
338 #endif
339
340 // Add an active element to the family tree.
341 if (elem->active())
342 family.push_back(elem);
343
344 // Or recurse into an ancestor element's children.
345 // Do not clear the vector any more.
346 else if (!elem->active())
347 for (auto & c : elem->child_ref_range())
348 if (!c.is_remote() && c.has_neighbor(neighbor_in))
349 ElemInternal::active_family_tree_by_neighbor (&c, family, neighbor_in, false);
350 }
351
352
353
354 template <class T>
355 void
356 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 libmesh_assert(!elem->subactive());
367 libmesh_assert(!elem->is_remote());
368
369 // Clear the vector if the flag reset tells us to.
370 if (reset)
371 family.clear();
372
373 // This only makes sense if we're already a topological neighbor
374 #ifndef NDEBUG
375 if (elem->level() >= neighbor_in->level())
376 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 if (elem->active())
384 family.push_back(elem);
385
386 // Or recurse into an ancestor element's children.
387 // Do not clear the vector any more.
388 else if (!elem->active())
389 for (auto & c : elem->child_ref_range())
390 if (!c.is_remote() &&
391 c.has_topological_neighbor(neighbor_in, mesh, point_locator, pb))
392 c.active_family_tree_by_topological_neighbor
393 (family, neighbor_in, mesh, point_locator, pb, false);
394 }
395
396
397
398 template <class T>
399 void
find_point_neighbors(T this_elem,std::set<T> & neighbor_set,T start_elem)400 find_point_neighbors(T this_elem,
401 std::set<T> & neighbor_set,
402 T start_elem)
403 {
404 libmesh_assert(start_elem);
405 libmesh_assert(start_elem->active());
406 libmesh_assert(start_elem->contains_vertex_of(this_elem) ||
407 this_elem->contains_vertex_of(start_elem));
408
409 neighbor_set.clear();
410 neighbor_set.insert(start_elem);
411
412 std::set<T> untested_set, next_untested_set;
413 untested_set.insert(start_elem);
414
415 while (!untested_set.empty())
416 {
417 // Loop over all the elements in the patch that haven't already
418 // been tested
419 for (const auto & elem : untested_set)
420 for (auto current_neighbor : elem->neighbor_ptr_range())
421 {
422 if (current_neighbor &&
423 !current_neighbor->is_remote()) // we have a real neighbor on this side
424 {
425 if (current_neighbor->active()) // ... if it is active
426 {
427 if (this_elem->contains_vertex_of(current_neighbor) // ... and touches us
428 || current_neighbor->contains_vertex_of(this_elem))
429 {
430 // Make sure we'll test it
431 if (!neighbor_set.count(current_neighbor))
432 next_untested_set.insert (current_neighbor);
433
434 // And add it
435 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 std::vector<T> active_neighbor_children;
443
444 current_neighbor->active_family_tree_by_neighbor
445 (active_neighbor_children, elem);
446
447 for (const auto & current_child : active_neighbor_children)
448 {
449 if (this_elem->contains_vertex_of(current_child) ||
450 current_child->contains_vertex_of(this_elem))
451 {
452 // Make sure we'll test it
453 if (!neighbor_set.count(current_child))
454 next_untested_set.insert (current_child);
455
456 neighbor_set.insert (current_child);
457 }
458 }
459 }
460 #endif // #ifdef LIBMESH_ENABLE_AMR
461 }
462 }
463 untested_set.swap(next_untested_set);
464 next_untested_set.clear();
465 }
466 }
467
468
469
470 template <class T>
471 void
find_interior_neighbors(T this_elem,std::set<T> & neighbor_set)472 find_interior_neighbors(T this_elem,
473 std::set<T> & neighbor_set)
474 {
475 neighbor_set.clear();
476
477 if ((this_elem->dim() >= LIBMESH_DIM) ||
478 !this_elem->interior_parent())
479 return;
480
481 T ip = this_elem->interior_parent();
482 libmesh_assert (ip->contains_vertex_of(this_elem) ||
483 this_elem->contains_vertex_of(ip));
484
485 libmesh_assert (!ip->subactive());
486
487 #ifdef LIBMESH_ENABLE_AMR
488 while (!ip->active()) // only possible with AMR, be careful because
489 { // ip->child_ptr(c) is only good with AMR.
490 for (auto & child : ip->child_ref_range())
491 {
492 if (child.contains_vertex_of(this_elem) ||
493 this_elem->contains_vertex_of(&child))
494 {
495 ip = &child;
496 break;
497 }
498 }
499 }
500 #endif
501
502 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 auto it = neighbor_set.begin();
514 const auto end = neighbor_set.end();
515
516 while (it != end)
517 {
518 auto current = it++;
519 T current_elem = *current;
520
521 // This won't invalidate iterator it, because it is already
522 // pointing to the next element
523 if (current_elem->level() > this_elem->level())
524 {
525 unsigned int vertices_contained = 0;
526 for (auto & n : current_elem->node_ref_range())
527 if (this_elem->contains_point(n))
528 vertices_contained++;
529
530 if (vertices_contained <= this_elem->dim())
531 {
532 neighbor_set.erase(current);
533 continue;
534 }
535 }
536 else
537 {
538 for (auto & n : this_elem->node_ref_range())
539 {
540 if (!current_elem->contains_point(n))
541 {
542 neighbor_set.erase(current);
543 break;
544 }
545 }
546 }
547 }
548 }
549
550 } // namespace ElemInternal
551 } // namespace libMesh
552
553 #endif
554