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
19
20 // C++ includes
21 #include <algorithm> // for std::sort
22 #include <array>
23 #include <iterator> // for std::ostream_iterator
24 #include <sstream>
25 #include <limits> // for std::numeric_limits<>
26 #include <cmath> // for std::sqrt()
27
28 // Local includes
29 #include "libmesh/auto_ptr.h" // libmesh_make_unique
30 #include "libmesh/elem.h"
31 #include "libmesh/fe_type.h"
32 #include "libmesh/fe_interface.h"
33 #include "libmesh/node_elem.h"
34 #include "libmesh/edge_edge2.h"
35 #include "libmesh/edge_edge3.h"
36 #include "libmesh/edge_edge4.h"
37 #include "libmesh/edge_inf_edge2.h"
38 #include "libmesh/face_tri3.h"
39 #include "libmesh/face_tri3_subdivision.h"
40 #include "libmesh/face_tri3_shell.h"
41 #include "libmesh/face_tri6.h"
42 #include "libmesh/face_quad4.h"
43 #include "libmesh/face_quad4_shell.h"
44 #include "libmesh/face_quad8.h"
45 #include "libmesh/face_quad8_shell.h"
46 #include "libmesh/face_quad9.h"
47 #include "libmesh/face_inf_quad4.h"
48 #include "libmesh/face_inf_quad6.h"
49 #include "libmesh/cell_tet4.h"
50 #include "libmesh/cell_tet10.h"
51 #include "libmesh/cell_hex8.h"
52 #include "libmesh/cell_hex20.h"
53 #include "libmesh/cell_hex27.h"
54 #include "libmesh/cell_inf_hex8.h"
55 #include "libmesh/cell_inf_hex16.h"
56 #include "libmesh/cell_inf_hex18.h"
57 #include "libmesh/cell_prism6.h"
58 #include "libmesh/cell_prism15.h"
59 #include "libmesh/cell_prism18.h"
60 #include "libmesh/cell_inf_prism6.h"
61 #include "libmesh/cell_inf_prism12.h"
62 #include "libmesh/cell_pyramid5.h"
63 #include "libmesh/cell_pyramid13.h"
64 #include "libmesh/cell_pyramid14.h"
65 #include "libmesh/fe_base.h"
66 #include "libmesh/mesh_base.h"
67 #include "libmesh/quadrature_gauss.h"
68 #include "libmesh/remote_elem.h"
69 #include "libmesh/reference_elem.h"
70 #include "libmesh/enum_to_string.h"
71 #include "libmesh/threads.h"
72 #include "libmesh/enum_elem_quality.h"
73 #include "libmesh/enum_io_package.h"
74 #include "libmesh/enum_order.h"
75 #include "libmesh/elem_internal.h"
76
77 #ifdef LIBMESH_ENABLE_PERIODIC
78 #include "libmesh/mesh.h"
79 #include "libmesh/periodic_boundaries.h"
80 #include "libmesh/boundary_info.h"
81 #endif
82
83 namespace libMesh
84 {
85
86 Threads::spin_mutex parent_indices_mutex;
87 Threads::spin_mutex parent_bracketing_nodes_mutex;
88
89 const subdomain_id_type Elem::invalid_subdomain_id = std::numeric_limits<subdomain_id_type>::max();
90
91 // Initialize static member variables
92 const unsigned int Elem::max_n_nodes;
93
94 const unsigned int Elem::type_to_n_nodes_map [] =
95 {
96 2, // EDGE2
97 3, // EDGE3
98 4, // EDGE4
99
100 3, // TRI3
101 6, // TRI6
102
103 4, // QUAD4
104 8, // QUAD8
105 9, // QUAD9
106
107 4, // TET4
108 10, // TET10
109
110 8, // HEX8
111 20, // HEX20
112 27, // HEX27
113
114 6, // PRISM6
115 15, // PRISM15
116 18, // PRISM18
117
118 5, // PYRAMID5
119 13, // PYRAMID13
120 14, // PYRAMID14
121
122 2, // INFEDGE2
123
124 4, // INFQUAD4
125 6, // INFQUAD6
126
127 8, // INFHEX8
128 16, // INFHEX16
129 18, // INFHEX18
130
131 6, // INFPRISM6
132 12, // INFPRISM12
133
134 1, // NODEELEM
135
136 0, // REMOTEELEM
137
138 3, // TRI3SUBDIVISION
139 3, // TRISHELL3
140 4, // QUADSHELL4
141 8, // QUADSHELL8
142 };
143
144 const unsigned int Elem::type_to_n_sides_map [] =
145 {
146 2, // EDGE2
147 2, // EDGE3
148 2, // EDGE4
149
150 3, // TRI3
151 3, // TRI6
152
153 4, // QUAD4
154 4, // QUAD8
155 4, // QUAD9
156
157 4, // TET4
158 4, // TET10
159
160 6, // HEX8
161 6, // HEX20
162 6, // HEX27
163
164 5, // PRISM6
165 5, // PRISM15
166 5, // PRISM18
167
168 5, // PYRAMID5
169 5, // PYRAMID13
170 5, // PYRAMID14
171
172 2, // INFEDGE2
173
174 3, // INFQUAD4
175 3, // INFQUAD6
176
177 5, // INFHEX8
178 5, // INFHEX16
179 5, // INFHEX18
180
181 4, // INFPRISM6
182 4, // INFPRISM12
183
184 0, // NODEELEM
185
186 0, // REMOTEELEM
187
188 3, // TRI3SUBDIVISION
189 3, // TRISHELL3
190 4, // QUADSHELL4
191 4, // QUADSHELL8
192 };
193
194 const unsigned int Elem::type_to_n_edges_map [] =
195 {
196 0, // EDGE2
197 0, // EDGE3
198 0, // EDGE4
199
200 3, // TRI3
201 3, // TRI6
202
203 4, // QUAD4
204 4, // QUAD8
205 4, // QUAD9
206
207 6, // TET4
208 6, // TET10
209
210 12, // HEX8
211 12, // HEX20
212 12, // HEX27
213
214 9, // PRISM6
215 9, // PRISM15
216 9, // PRISM18
217
218 8, // PYRAMID5
219 8, // PYRAMID13
220 8, // PYRAMID14
221
222 0, // INFEDGE2
223
224 4, // INFQUAD4
225 4, // INFQUAD6
226
227 8, // INFHEX8
228 8, // INFHEX16
229 8, // INFHEX18
230
231 6, // INFPRISM6
232 6, // INFPRISM12
233
234 0, // NODEELEM
235
236 0, // REMOTEELEM
237
238 3, // TRI3SUBDIVISION
239 3, // TRISHELL3
240 4, // QUADSHELL4
241 4, // QUADSHELL8
242 };
243
244 // ------------------------------------------------------------
245 // Elem class member functions
build(const ElemType type,Elem * p)246 std::unique_ptr<Elem> Elem::build(const ElemType type,
247 Elem * p)
248 {
249 switch (type)
250 {
251 // 0D elements
252 case NODEELEM:
253 return libmesh_make_unique<NodeElem>(p);
254
255 // 1D elements
256 case EDGE2:
257 return libmesh_make_unique<Edge2>(p);
258 case EDGE3:
259 return libmesh_make_unique<Edge3>(p);
260 case EDGE4:
261 return libmesh_make_unique<Edge4>(p);
262
263 // 2D elements
264 case TRI3:
265 return libmesh_make_unique<Tri3>(p);
266 case TRISHELL3:
267 return libmesh_make_unique<TriShell3>(p);
268 case TRI3SUBDIVISION:
269 return libmesh_make_unique<Tri3Subdivision>(p);
270 case TRI6:
271 return libmesh_make_unique<Tri6>(p);
272 case QUAD4:
273 return libmesh_make_unique<Quad4>(p);
274 case QUADSHELL4:
275 return libmesh_make_unique<QuadShell4>(p);
276 case QUAD8:
277 return libmesh_make_unique<Quad8>(p);
278 case QUADSHELL8:
279 return libmesh_make_unique<QuadShell8>(p);
280 case QUAD9:
281 return libmesh_make_unique<Quad9>(p);
282
283 // 3D elements
284 case TET4:
285 return libmesh_make_unique<Tet4>(p);
286 case TET10:
287 return libmesh_make_unique<Tet10>(p);
288 case HEX8:
289 return libmesh_make_unique<Hex8>(p);
290 case HEX20:
291 return libmesh_make_unique<Hex20>(p);
292 case HEX27:
293 return libmesh_make_unique<Hex27>(p);
294 case PRISM6:
295 return libmesh_make_unique<Prism6>(p);
296 case PRISM15:
297 return libmesh_make_unique<Prism15>(p);
298 case PRISM18:
299 return libmesh_make_unique<Prism18>(p);
300 case PYRAMID5:
301 return libmesh_make_unique<Pyramid5>(p);
302 case PYRAMID13:
303 return libmesh_make_unique<Pyramid13>(p);
304 case PYRAMID14:
305 return libmesh_make_unique<Pyramid14>(p);
306
307 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
308 // 1D infinite elements
309 case INFEDGE2:
310 return libmesh_make_unique<InfEdge2>(p);
311
312 // 2D infinite elements
313 case INFQUAD4:
314 return libmesh_make_unique<InfQuad4>(p);
315 case INFQUAD6:
316 return libmesh_make_unique<InfQuad6>(p);
317
318 // 3D infinite elements
319 case INFHEX8:
320 return libmesh_make_unique<InfHex8>(p);
321 case INFHEX16:
322 return libmesh_make_unique<InfHex16>(p);
323 case INFHEX18:
324 return libmesh_make_unique<InfHex18>(p);
325 case INFPRISM6:
326 return libmesh_make_unique<InfPrism6>(p);
327 case INFPRISM12:
328 return libmesh_make_unique<InfPrism12>(p);
329 #endif
330
331 default:
332 libmesh_error_msg("ERROR: Undefined element type!");
333 }
334 }
335
336
337
build_with_id(const ElemType type,dof_id_type id)338 std::unique_ptr<Elem> Elem::build_with_id (const ElemType type,
339 dof_id_type id)
340 {
341 // Call the other build() method with nullptr parent, then set the
342 // required id.
343 auto temp = Elem::build(type, nullptr);
344 temp->set_id(id);
345 return temp;
346 }
347
348
349
reference_elem()350 const Elem * Elem::reference_elem () const
351 {
352 return &(ReferenceElem::get(this->type()));
353 }
354
355
356
centroid()357 Point Elem::centroid() const
358 {
359 Point cp;
360
361 const auto n_vertices = this->n_vertices();
362
363 for (unsigned int n=0; n<n_vertices; n++)
364 cp.add (this->point(n));
365
366 return (cp /= static_cast<Real>(n_vertices));
367 }
368
369
370
hmin()371 Real Elem::hmin() const
372 {
373 Real h_min=std::numeric_limits<Real>::max();
374
375 // Avoid calling a virtual a lot of times
376 const auto n_vertices = this->n_vertices();
377
378 for (unsigned int n_outer=0; n_outer<n_vertices; n_outer++)
379 for (unsigned int n_inner=n_outer+1; n_inner<n_vertices; n_inner++)
380 {
381 const auto diff = (this->point(n_outer) - this->point(n_inner));
382
383 h_min = std::min(h_min, diff.norm_sq());
384 }
385
386 return std::sqrt(h_min);
387 }
388
389
390
hmax()391 Real Elem::hmax() const
392 {
393 Real h_max=0;
394
395 // Avoid calling a virtual a lot of times
396 const auto n_vertices = this->n_vertices();
397
398 for (unsigned int n_outer=0; n_outer<n_vertices; n_outer++)
399 for (unsigned int n_inner=n_outer+1; n_inner<n_vertices; n_inner++)
400 {
401 const auto diff = (this->point(n_outer) - this->point(n_inner));
402
403 h_max = std::max(h_max, diff.norm_sq());
404 }
405
406 return std::sqrt(h_max);
407 }
408
409
410
length(const unsigned int n1,const unsigned int n2)411 Real Elem::length(const unsigned int n1,
412 const unsigned int n2) const
413 {
414 libmesh_assert_less ( n1, this->n_vertices() );
415 libmesh_assert_less ( n2, this->n_vertices() );
416
417 return (this->point(n1) - this->point(n2)).norm();
418 }
419
420
421
key()422 dof_id_type Elem::key () const
423 {
424 const unsigned short n_n = this->n_nodes();
425
426 std::array<dof_id_type, Elem::max_n_nodes> node_ids;
427
428 for (unsigned short n=0; n != n_n; ++n)
429 node_ids[n] = this->node_id(n);
430
431 // Always sort, so that different local node numberings hash to the
432 // same value.
433 std::sort (node_ids.begin(), node_ids.begin()+n_n);
434
435 return Utility::hashword(node_ids.data(), n_n);
436 }
437
438
439
440 bool Elem::operator == (const Elem & rhs) const
441 {
442 // If the elements aren't the same type, they aren't equal
443 if (this->type() != rhs.type())
444 return false;
445
446 const unsigned short n_n = this->n_nodes();
447 libmesh_assert_equal_to(n_n, rhs.n_nodes());
448
449 // Make two sorted arrays of global node ids and compare them for
450 // equality.
451 std::array<dof_id_type, Elem::max_n_nodes> this_ids, rhs_ids;
452
453 for (unsigned short n = 0; n != n_n; n++)
454 {
455 this_ids[n] = this->node_id(n);
456 rhs_ids[n] = rhs.node_id(n);
457 }
458
459 // Sort the vectors to rule out different local node numberings.
460 std::sort(this_ids.begin(), this_ids.begin()+n_n);
461 std::sort(rhs_ids.begin(), rhs_ids.begin()+n_n);
462
463 // If the node ids match, the elements are equal!
464 for (unsigned short n = 0; n != n_n; ++n)
465 if (this_ids[n] != rhs_ids[n])
466 return false;
467 return true;
468 }
469
470
471
is_semilocal(const processor_id_type my_pid)472 bool Elem::is_semilocal(const processor_id_type my_pid) const
473 {
474 std::set<const Elem *> point_neighbors;
475
476 this->find_point_neighbors(point_neighbors);
477
478 for (const auto & elem : point_neighbors)
479 if (elem->processor_id() == my_pid)
480 return true;
481
482 return false;
483 }
484
485
486
which_side_am_i(const Elem * e)487 unsigned int Elem::which_side_am_i (const Elem * e) const
488 {
489 libmesh_assert(e);
490
491 const unsigned int ns = this->n_sides();
492 const unsigned int nn = this->n_nodes();
493
494 const unsigned int en = e->n_nodes();
495
496 // e might be on any side until proven otherwise
497 std::vector<bool> might_be_side(ns, true);
498
499 for (unsigned int i=0; i != en; ++i)
500 {
501 Point side_point = e->point(i);
502 unsigned int local_node_id = libMesh::invalid_uint;
503
504 // Look for a node of this that's contiguous with node i of
505 // e. Note that the exact floating point comparison of Point
506 // positions is intentional, see the class documentation for
507 // this function.
508 for (unsigned int j=0; j != nn; ++j)
509 if (this->point(j) == side_point)
510 local_node_id = j;
511
512 // If a node of e isn't contiguous with some node of this, then
513 // e isn't a side of this.
514 if (local_node_id == libMesh::invalid_uint)
515 return libMesh::invalid_uint;
516
517 // If a node of e isn't contiguous with some node on side s of
518 // this, then e isn't on side s.
519 for (unsigned int s=0; s != ns; ++s)
520 if (!this->is_node_on_side(local_node_id, s))
521 might_be_side[s] = false;
522 }
523
524 for (unsigned int s=0; s != ns; ++s)
525 if (might_be_side[s])
526 {
527 #ifdef DEBUG
528 for (unsigned int s2=s+1; s2 < ns; ++s2)
529 libmesh_assert (!might_be_side[s2]);
530 #endif
531 return s;
532 }
533
534 // Didn't find any matching side
535 return libMesh::invalid_uint;
536 }
537
538
539
which_node_am_i(unsigned int side,unsigned int side_node)540 unsigned int Elem::which_node_am_i(unsigned int side,
541 unsigned int side_node) const
542 {
543 libmesh_deprecated();
544 return local_side_node(side, side_node);
545 }
546
547
548
contains_vertex_of(const Elem * e)549 bool Elem::contains_vertex_of(const Elem * e) const
550 {
551 // Our vertices are the first numbered nodes
552 for (auto n : make_range(e->n_vertices()))
553 if (this->contains_point(e->point(n)))
554 return true;
555 return false;
556 }
557
558
559
contains_edge_of(const Elem * e)560 bool Elem::contains_edge_of(const Elem * e) const
561 {
562 unsigned int num_contained_edges = 0;
563
564 // Our vertices are the first numbered nodes
565 for (auto n : make_range(e->n_vertices()))
566 {
567 if (this->contains_point(e->point(n)))
568 {
569 num_contained_edges++;
570 if (num_contained_edges>=2)
571 {
572 return true;
573 }
574 }
575 }
576 return false;
577 }
578
579
580
find_point_neighbors(const Point & p,std::set<const Elem * > & neighbor_set)581 void Elem::find_point_neighbors(const Point & p,
582 std::set<const Elem *> & neighbor_set) const
583 {
584 libmesh_assert(this->contains_point(p));
585 libmesh_assert(this->active());
586
587 neighbor_set.clear();
588 neighbor_set.insert(this);
589
590 std::set<const Elem *> untested_set, next_untested_set;
591 untested_set.insert(this);
592
593 #ifdef LIBMESH_ENABLE_AMR
594 std::vector<const Elem *> active_neighbor_children;
595 #endif // #ifdef LIBMESH_ENABLE_AMR
596
597 while (!untested_set.empty())
598 {
599 // Loop over all the elements in the patch that haven't already
600 // been tested
601 for (const auto & elem : untested_set)
602 for (auto current_neighbor : elem->neighbor_ptr_range())
603 {
604 if (current_neighbor &&
605 current_neighbor != remote_elem) // we have a real neighbor on this side
606 {
607 if (current_neighbor->active()) // ... if it is active
608 {
609 auto it = neighbor_set.lower_bound(current_neighbor);
610 if ((it == neighbor_set.end() || *it != current_neighbor) &&
611 current_neighbor->contains_point(p)) // ... don't have and touches p
612 {
613 // Add it and test it
614 next_untested_set.insert(current_neighbor);
615 neighbor_set.emplace_hint(it, current_neighbor);
616 }
617 }
618 #ifdef LIBMESH_ENABLE_AMR
619 else // ... the neighbor is *not* active,
620 { // ... so add *all* neighboring
621 // active children that touch p
622 active_neighbor_children.clear();
623 current_neighbor->active_family_tree_by_neighbor
624 (active_neighbor_children, elem);
625
626 for (const auto & current_child : active_neighbor_children)
627 {
628 auto it = neighbor_set.lower_bound(current_child);
629 if ((it == neighbor_set.end() || *it != current_child) &&
630 current_child->contains_point(p))
631 {
632 // Add it and test it
633 next_untested_set.insert(current_child);
634 neighbor_set.emplace_hint(it, current_child);
635 }
636 }
637 }
638 #endif // #ifdef LIBMESH_ENABLE_AMR
639 }
640 }
641 untested_set.swap(next_untested_set);
642 next_untested_set.clear();
643 }
644 }
645
646
647
find_point_neighbors(std::set<const Elem * > & neighbor_set)648 void Elem::find_point_neighbors(std::set<const Elem *> & neighbor_set) const
649 {
650 this->find_point_neighbors(neighbor_set, this);
651 }
652
653
654
find_point_neighbors(std::set<const Elem * > & neighbor_set,const Elem * start_elem)655 void Elem::find_point_neighbors(std::set<const Elem *> & neighbor_set,
656 const Elem * start_elem) const
657 {
658 ElemInternal::find_point_neighbors(this, neighbor_set, start_elem);
659 }
660
661
662
find_point_neighbors(std::set<Elem * > & neighbor_set,Elem * start_elem)663 void Elem::find_point_neighbors(std::set<Elem *> & neighbor_set,
664 Elem * start_elem)
665 {
666 ElemInternal::find_point_neighbors(this, neighbor_set, start_elem);
667 }
668
669
670
find_edge_neighbors(const Point & p1,const Point & p2,std::set<const Elem * > & neighbor_set)671 void Elem::find_edge_neighbors(const Point & p1,
672 const Point & p2,
673 std::set<const Elem *> & neighbor_set) const
674 {
675 // Simple but perhaps suboptimal code: find elements containing the
676 // first point, then winnow this set down by removing elements which
677 // don't also contain the second point
678
679 libmesh_assert(this->contains_point(p2));
680 this->find_point_neighbors(p1, neighbor_set);
681
682 std::set<const Elem *>::iterator it = neighbor_set.begin();
683 const std::set<const Elem *>::iterator end = neighbor_set.end();
684
685 while (it != end)
686 {
687 // As of C++11, set::erase returns an iterator to the element
688 // following the erased element, or end.
689 if (!(*it)->contains_point(p2))
690 it = neighbor_set.erase(it);
691 else
692 ++it;
693 }
694 }
695
696
697
find_edge_neighbors(std::set<const Elem * > & neighbor_set)698 void Elem::find_edge_neighbors(std::set<const Elem *> & neighbor_set) const
699 {
700 neighbor_set.clear();
701 neighbor_set.insert(this);
702
703 std::set<const Elem *> untested_set, next_untested_set;
704 untested_set.insert(this);
705
706 while (!untested_set.empty())
707 {
708 // Loop over all the elements in the patch that haven't already
709 // been tested
710 for (const auto & elem : untested_set)
711 {
712 for (auto current_neighbor : elem->neighbor_ptr_range())
713 {
714 if (current_neighbor &&
715 current_neighbor != remote_elem) // we have a real neighbor on this side
716 {
717 if (current_neighbor->active()) // ... if it is active
718 {
719 if (this->contains_edge_of(current_neighbor) // ... and touches us
720 || current_neighbor->contains_edge_of(this))
721 {
722 // Make sure we'll test it
723 if (!neighbor_set.count(current_neighbor))
724 next_untested_set.insert (current_neighbor);
725
726 // And add it
727 neighbor_set.insert (current_neighbor);
728 }
729 }
730 #ifdef LIBMESH_ENABLE_AMR
731 else // ... the neighbor is *not* active,
732 { // ... so add *all* neighboring
733 // active children
734 std::vector<const Elem *> active_neighbor_children;
735
736 current_neighbor->active_family_tree_by_neighbor
737 (active_neighbor_children, elem);
738
739 for (const auto & current_child : active_neighbor_children)
740 if (this->contains_edge_of(current_child) || current_child->contains_edge_of(this))
741 {
742 // Make sure we'll test it
743 if (!neighbor_set.count(current_child))
744 next_untested_set.insert (current_child);
745
746 neighbor_set.insert (current_child);
747 }
748 }
749 #endif // #ifdef LIBMESH_ENABLE_AMR
750 }
751 }
752 }
753 untested_set.swap(next_untested_set);
754 next_untested_set.clear();
755 }
756 }
757
758
759
find_interior_neighbors(std::set<const Elem * > & neighbor_set)760 void Elem::find_interior_neighbors(std::set<const Elem *> & neighbor_set) const
761 {
762 ElemInternal::find_interior_neighbors(this, neighbor_set);
763 }
764
765
766
find_interior_neighbors(std::set<Elem * > & neighbor_set)767 void Elem::find_interior_neighbors(std::set<Elem *> & neighbor_set)
768 {
769 ElemInternal::find_interior_neighbors(this, neighbor_set);
770 }
771
772
773
interior_parent()774 const Elem * Elem::interior_parent () const
775 {
776 // interior parents make no sense for full-dimensional elements.
777 if (this->dim() >= LIBMESH_DIM)
778 return nullptr;
779
780 // they USED TO BE only good for level-0 elements, but we now
781 // support keeping interior_parent() valid on refined boundary
782 // elements.
783 // if (this->level() != 0)
784 // return this->parent()->interior_parent();
785
786 // We store the interior_parent pointer after both the parent
787 // neighbor and neighbor pointers
788 Elem * interior_p = _elemlinks[1+this->n_sides()];
789
790 // If we have an interior_parent, we USED TO assume it was a
791 // one-higher-dimensional interior element, but we now allow e.g.
792 // edge elements to have a 3D interior_parent with no
793 // intermediate 2D element.
794 // libmesh_assert (!interior_p ||
795 // interior_p->dim() == (this->dim()+1));
796 libmesh_assert (!interior_p ||
797 (interior_p == remote_elem) ||
798 (interior_p->dim() > this->dim()));
799
800 return interior_p;
801 }
802
803
804
interior_parent()805 Elem * Elem::interior_parent ()
806 {
807 // See the const version for comments
808 if (this->dim() >= LIBMESH_DIM)
809 return nullptr;
810
811 Elem * interior_p = _elemlinks[1+this->n_sides()];
812
813 libmesh_assert (!interior_p ||
814 (interior_p == remote_elem) ||
815 (interior_p->dim() > this->dim()));
816
817 return interior_p;
818 }
819
820
821
set_interior_parent(Elem * p)822 void Elem::set_interior_parent (Elem * p)
823 {
824 // interior parents make no sense for full-dimensional elements.
825 libmesh_assert_less (this->dim(), LIBMESH_DIM);
826
827 // If we have an interior_parent, we USED TO assume it was a
828 // one-higher-dimensional interior element, but we now allow e.g.
829 // edge elements to have a 3D interior_parent with no
830 // intermediate 2D element.
831 // libmesh_assert (!p ||
832 // p->dim() == (this->dim()+1));
833 libmesh_assert (!p ||
834 (p == remote_elem) ||
835 (p->dim() > this->dim()));
836
837 _elemlinks[1+this->n_sides()] = p;
838 }
839
840
841
842 #ifdef LIBMESH_ENABLE_PERIODIC
843
topological_neighbor(const unsigned int i,MeshBase & mesh,const PointLocatorBase & point_locator,const PeriodicBoundaries * pb)844 Elem * Elem::topological_neighbor (const unsigned int i,
845 MeshBase & mesh,
846 const PointLocatorBase & point_locator,
847 const PeriodicBoundaries * pb)
848 {
849 libmesh_assert_less (i, this->n_neighbors());
850
851 Elem * neighbor_i = this->neighbor_ptr(i);
852 if (neighbor_i != nullptr)
853 return neighbor_i;
854
855 if (pb)
856 {
857 // Since the neighbor is nullptr it must be on a boundary. We need
858 // see if this is a periodic boundary in which case it will have a
859 // topological neighbor
860 std::vector<boundary_id_type> bc_ids;
861 mesh.get_boundary_info().boundary_ids(this, cast_int<unsigned short>(i), bc_ids);
862 for (const auto & id : bc_ids)
863 if (pb->boundary(id))
864 {
865 // Since the point locator inside of periodic boundaries
866 // returns a const pointer we will retrieve the proper
867 // pointer directly from the mesh object.
868 const Elem * const cn = pb->neighbor(id, point_locator, this, i);
869 neighbor_i = const_cast<Elem *>(cn);
870
871 // Since coarse elements do not have more refined
872 // neighbors we need to make sure that we don't return one
873 // of these types of neighbors.
874 if (neighbor_i)
875 while (level() < neighbor_i->level())
876 neighbor_i = neighbor_i->parent();
877 return neighbor_i;
878 }
879 }
880
881 return nullptr;
882 }
883
884
885
topological_neighbor(const unsigned int i,const MeshBase & mesh,const PointLocatorBase & point_locator,const PeriodicBoundaries * pb)886 const Elem * Elem::topological_neighbor (const unsigned int i,
887 const MeshBase & mesh,
888 const PointLocatorBase & point_locator,
889 const PeriodicBoundaries * pb) const
890 {
891 libmesh_assert_less (i, this->n_neighbors());
892
893 const Elem * neighbor_i = this->neighbor_ptr(i);
894 if (neighbor_i != nullptr)
895 return neighbor_i;
896
897 if (pb)
898 {
899 // Since the neighbor is nullptr it must be on a boundary. We need
900 // see if this is a periodic boundary in which case it will have a
901 // topological neighbor
902 std::vector<boundary_id_type> bc_ids;
903 mesh.get_boundary_info().boundary_ids(this, cast_int<unsigned short>(i), bc_ids);
904 for (const auto & id : bc_ids)
905 if (pb->boundary(id))
906 {
907 neighbor_i = pb->neighbor(id, point_locator, this, i);
908
909 // Since coarse elements do not have more refined
910 // neighbors we need to make sure that we don't return one
911 // of these types of neighbors.
912 if (neighbor_i)
913 while (level() < neighbor_i->level())
914 neighbor_i = neighbor_i->parent();
915 return neighbor_i;
916 }
917 }
918
919 return nullptr;
920 }
921
922
has_topological_neighbor(const Elem * elem,const MeshBase & mesh,const PointLocatorBase & point_locator,const PeriodicBoundaries * pb)923 bool Elem::has_topological_neighbor (const Elem * elem,
924 const MeshBase & mesh,
925 const PointLocatorBase & point_locator,
926 const PeriodicBoundaries * pb) const
927 {
928 // First see if this is a normal "interior" neighbor
929 if (has_neighbor(elem))
930 return true;
931
932 for (auto n : this->side_index_range())
933 if (this->topological_neighbor(n, mesh, point_locator, pb))
934 return true;
935
936 return false;
937 }
938
939
940 #endif
941
942 #ifdef DEBUG
943
libmesh_assert_valid_node_pointers()944 void Elem::libmesh_assert_valid_node_pointers() const
945 {
946 libmesh_assert(this->valid_id());
947 for (auto n : this->node_index_range())
948 {
949 libmesh_assert(this->node_ptr(n));
950 libmesh_assert(this->node_ptr(n)->valid_id());
951 }
952 }
953
954
955
libmesh_assert_valid_neighbors()956 void Elem::libmesh_assert_valid_neighbors() const
957 {
958 for (auto n : this->side_index_range())
959 {
960 const Elem * neigh = this->neighbor_ptr(n);
961
962 // Any element might have a remote neighbor; checking
963 // to make sure that's not inaccurate is tough.
964 if (neigh == remote_elem)
965 continue;
966
967 if (neigh)
968 {
969 // Only subactive elements have subactive neighbors
970 libmesh_assert (this->subactive() || !neigh->subactive());
971
972 const Elem * elem = this;
973
974 // If we're subactive but our neighbor isn't, its
975 // return neighbor link will be to our first active
976 // ancestor OR to our inactive ancestor of the same
977 // level as neigh,
978 if (this->subactive() && !neigh->subactive())
979 {
980 for (elem = this; !elem->active();
981 elem = elem->parent())
982 libmesh_assert(elem);
983 }
984 else
985 {
986 unsigned int rev = neigh->which_neighbor_am_i(elem);
987 libmesh_assert_less (rev, neigh->n_neighbors());
988
989 if (this->subactive() && !neigh->subactive())
990 {
991 while (neigh->neighbor_ptr(rev) != elem)
992 {
993 libmesh_assert(elem->parent());
994 elem = elem->parent();
995 }
996 }
997 else
998 {
999 const Elem * nn = neigh->neighbor_ptr(rev);
1000 libmesh_assert(nn);
1001
1002 for (; elem != nn; elem = elem->parent())
1003 libmesh_assert(elem);
1004 }
1005 }
1006 }
1007 // If we don't have a neighbor and we're not subactive, our
1008 // ancestors shouldn't have any neighbors in this same
1009 // direction.
1010 else if (!this->subactive())
1011 {
1012 const Elem * my_parent = this->parent();
1013 if (my_parent &&
1014 // A parent with a different dimension isn't really one of
1015 // our ancestors, it means we're on a boundary mesh and this
1016 // is an interior mesh element for which we're on a side.
1017 // Nothing to test for in that case.
1018 (my_parent->dim() == this->dim()))
1019 libmesh_assert (!my_parent->neighbor_ptr(n));
1020 }
1021 }
1022 }
1023
1024 #endif // DEBUG
1025
1026
1027
make_links_to_me_local(unsigned int n,unsigned int nn)1028 void Elem::make_links_to_me_local(unsigned int n, unsigned int nn)
1029 {
1030 Elem * neigh = this->neighbor_ptr(n);
1031
1032 // Don't bother calling this function unless it's necessary
1033 libmesh_assert(neigh);
1034 libmesh_assert(!neigh->is_remote());
1035
1036 // We never have neighbors more refined than us
1037 libmesh_assert_less_equal (neigh->level(), this->level());
1038
1039 // We never have subactive neighbors of non subactive elements
1040 libmesh_assert(!neigh->subactive() || this->subactive());
1041
1042 // If we have a neighbor less refined than us then it must not
1043 // have any more refined descendants we could have pointed to
1044 // instead.
1045 libmesh_assert((neigh->level() == this->level()) ||
1046 (neigh->active() && !this->subactive()) ||
1047 (!neigh->has_children() && this->subactive()));
1048
1049 // If neigh is at our level, then its family might have
1050 // remote_elem neighbor links which need to point to us
1051 // instead, but if not, then we're done.
1052 if (neigh->level() != this->level())
1053 return;
1054
1055 // What side of neigh are we on? nn.
1056 //
1057 // We can't use the usual Elem method because we're in the middle of
1058 // restoring topology. We can't compare side_ptr nodes because
1059 // users want to abuse neighbor_ptr to point to
1060 // not-technically-neighbors across mesh slits. We can't compare
1061 // node locations because users want to move those
1062 // not-technically-neighbors until they're
1063 // not-even-geometrically-neighbors.
1064
1065 // Find any elements that ought to point to elem
1066 std::vector<Elem *> neigh_family;
1067 #ifdef LIBMESH_ENABLE_AMR
1068 if (this->active())
1069 neigh->family_tree_by_side(neigh_family, nn);
1070 else
1071 #endif
1072 neigh_family.push_back(neigh);
1073
1074 // And point them to elem
1075 for (auto & neigh_family_member : neigh_family)
1076 {
1077 // Only subactive elements point to other subactive elements
1078 if (this->subactive() && !neigh_family_member->subactive())
1079 continue;
1080
1081 // Ideally, the neighbor link ought to either be correct
1082 // already or ought to be to remote_elem.
1083 //
1084 // However, if we're redistributing a newly created elem,
1085 // after an AMR step but before find_neighbors has fixed up
1086 // neighbor links, we might have an out of date neighbor
1087 // link to elem's parent instead.
1088 #ifdef LIBMESH_ENABLE_AMR
1089 libmesh_assert((neigh_family_member->neighbor_ptr(nn) &&
1090 (neigh_family_member->neighbor_ptr(nn)->active() ||
1091 neigh_family_member->neighbor_ptr(nn)->is_ancestor_of(this))) ||
1092 (neigh_family_member->neighbor_ptr(nn) == remote_elem) ||
1093 ((this->refinement_flag() == JUST_REFINED) &&
1094 (this->parent() != nullptr) &&
1095 (neigh_family_member->neighbor_ptr(nn) == this->parent())));
1096 #else
1097 libmesh_assert((neigh_family_member->neighbor_ptr(nn) == this) ||
1098 (neigh_family_member->neighbor_ptr(nn) == remote_elem));
1099 #endif
1100
1101 neigh_family_member->set_neighbor(nn, this);
1102 }
1103 }
1104
1105
make_links_to_me_remote()1106 void Elem::make_links_to_me_remote()
1107 {
1108 libmesh_assert_not_equal_to (this, remote_elem);
1109
1110 // We need to have handled any children first
1111 #if defined(LIBMESH_ENABLE_AMR) && defined(DEBUG)
1112 if (this->has_children())
1113 for (auto & child : this->child_ref_range())
1114 libmesh_assert_equal_to (&child, remote_elem);
1115 #endif
1116
1117 // Remotify any neighbor links
1118 for (auto neigh : this->neighbor_ptr_range())
1119 {
1120 if (neigh && neigh != remote_elem)
1121 {
1122 // My neighbor should never be more refined than me; my real
1123 // neighbor would have been its parent in that case.
1124 libmesh_assert_greater_equal (this->level(), neigh->level());
1125
1126 if (this->level() == neigh->level() &&
1127 neigh->has_neighbor(this))
1128 {
1129 #ifdef LIBMESH_ENABLE_AMR
1130 // My neighbor may have descendants which also consider me a
1131 // neighbor
1132 std::vector<Elem *> family;
1133 neigh->total_family_tree_by_neighbor (family, this);
1134
1135 // FIXME - There's a lot of ugly const_casts here; we
1136 // may want to make remote_elem non-const
1137 for (auto & n : family)
1138 {
1139 libmesh_assert (n);
1140 if (n == remote_elem)
1141 continue;
1142 unsigned int my_s = n->which_neighbor_am_i(this);
1143 libmesh_assert_less (my_s, n->n_neighbors());
1144 libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1145 n->set_neighbor(my_s, const_cast<RemoteElem *>(remote_elem));
1146 }
1147 #else
1148 unsigned int my_s = neigh->which_neighbor_am_i(this);
1149 libmesh_assert_less (my_s, neigh->n_neighbors());
1150 libmesh_assert_equal_to (neigh->neighbor_ptr(my_s), this);
1151 neigh->set_neighbor(my_s, const_cast<RemoteElem *>(remote_elem));
1152 #endif
1153 }
1154 #ifdef LIBMESH_ENABLE_AMR
1155 // Even if my neighbor doesn't link back to me, it might
1156 // have subactive descendants which do
1157 else if (neigh->has_children())
1158 {
1159 // If my neighbor at the same level doesn't have me as a
1160 // neighbor, I must be subactive
1161 libmesh_assert(this->level() > neigh->level() ||
1162 this->subactive());
1163
1164 // My neighbor must have some ancestor of mine as a
1165 // neighbor
1166 Elem * my_ancestor = this->parent();
1167 libmesh_assert(my_ancestor);
1168 while (!neigh->has_neighbor(my_ancestor))
1169 {
1170 my_ancestor = my_ancestor->parent();
1171 libmesh_assert(my_ancestor);
1172 }
1173
1174 // My neighbor may have descendants which consider me a
1175 // neighbor
1176 std::vector<Elem *> family;
1177 neigh->total_family_tree_by_subneighbor (family, my_ancestor, this);
1178
1179 for (auto & n : family)
1180 {
1181 libmesh_assert (n);
1182 if (n->is_remote())
1183 continue;
1184 unsigned int my_s = n->which_neighbor_am_i(this);
1185 libmesh_assert_less (my_s, n->n_neighbors());
1186 libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1187 // TODO: we may want to make remote_elem non-const.
1188 n->set_neighbor(my_s, const_cast<RemoteElem *>(remote_elem));
1189 }
1190 }
1191 #endif
1192 }
1193 }
1194
1195 #ifdef LIBMESH_ENABLE_AMR
1196 // Remotify parent's child link
1197 Elem * my_parent = this->parent();
1198 if (my_parent &&
1199 // As long as it's not already remote
1200 my_parent != remote_elem &&
1201 // And it's a real parent, not an interior parent
1202 this->dim() == my_parent->dim())
1203 {
1204 unsigned int me = my_parent->which_child_am_i(this);
1205 libmesh_assert_equal_to (my_parent->child_ptr(me), this);
1206 my_parent->set_child(me, const_cast<RemoteElem *>(remote_elem));
1207 }
1208 #endif
1209 }
1210
1211
remove_links_to_me()1212 void Elem::remove_links_to_me()
1213 {
1214 libmesh_assert_not_equal_to (this, remote_elem);
1215
1216 // We need to have handled any children first
1217 #ifdef LIBMESH_ENABLE_AMR
1218 libmesh_assert (!this->has_children());
1219 #endif
1220
1221 // Nullify any neighbor links
1222 for (auto neigh : this->neighbor_ptr_range())
1223 {
1224 if (neigh && neigh != remote_elem)
1225 {
1226 // My neighbor should never be more refined than me; my real
1227 // neighbor would have been its parent in that case.
1228 libmesh_assert_greater_equal (this->level(), neigh->level());
1229
1230 if (this->level() == neigh->level() &&
1231 neigh->has_neighbor(this))
1232 {
1233 #ifdef LIBMESH_ENABLE_AMR
1234 // My neighbor may have descendants which also consider me a
1235 // neighbor
1236 std::vector<Elem *> family;
1237 neigh->total_family_tree_by_neighbor (family, this);
1238
1239 for (auto & n : family)
1240 {
1241 libmesh_assert (n);
1242 if (n->is_remote())
1243 continue;
1244 unsigned int my_s = n->which_neighbor_am_i(this);
1245 libmesh_assert_less (my_s, n->n_neighbors());
1246 libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1247 n->set_neighbor(my_s, nullptr);
1248 }
1249 #else
1250 unsigned int my_s = neigh->which_neighbor_am_i(this);
1251 libmesh_assert_less (my_s, neigh->n_neighbors());
1252 libmesh_assert_equal_to (neigh->neighbor_ptr(my_s), this);
1253 neigh->set_neighbor(my_s, nullptr);
1254 #endif
1255 }
1256 #ifdef LIBMESH_ENABLE_AMR
1257 // Even if my neighbor doesn't link back to me, it might
1258 // have subactive descendants which do
1259 else if (neigh->has_children())
1260 {
1261 // If my neighbor at the same level doesn't have me as a
1262 // neighbor, I must be subactive
1263 libmesh_assert(this->level() > neigh->level() ||
1264 this->subactive());
1265
1266 // My neighbor must have some ancestor of mine as a
1267 // neighbor
1268 Elem * my_ancestor = this->parent();
1269 libmesh_assert(my_ancestor);
1270 while (!neigh->has_neighbor(my_ancestor))
1271 {
1272 my_ancestor = my_ancestor->parent();
1273 libmesh_assert(my_ancestor);
1274 }
1275
1276 // My neighbor may have descendants which consider me a
1277 // neighbor
1278 std::vector<Elem *> family;
1279 neigh->total_family_tree_by_subneighbor (family, my_ancestor, this);
1280
1281 for (auto & n : family)
1282 {
1283 libmesh_assert (n);
1284 if (n->is_remote())
1285 continue;
1286 unsigned int my_s = n->which_neighbor_am_i(this);
1287 libmesh_assert_less (my_s, n->n_neighbors());
1288 libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1289 n->set_neighbor(my_s, nullptr);
1290 }
1291 }
1292 #endif
1293 }
1294 }
1295
1296 #ifdef LIBMESH_ENABLE_AMR
1297 // We can't currently delete a child with a parent!
1298 libmesh_assert (!this->parent());
1299 #endif
1300 }
1301
1302
1303
write_connectivity(std::ostream & out_stream,const IOPackage iop)1304 void Elem::write_connectivity (std::ostream & out_stream,
1305 const IOPackage iop) const
1306 {
1307 libmesh_assert (out_stream.good());
1308 libmesh_assert(_nodes);
1309 libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
1310
1311 switch (iop)
1312 {
1313 case TECPLOT:
1314 {
1315 // This connectivity vector will be used repeatedly instead
1316 // of being reconstructed inside the loop.
1317 std::vector<dof_id_type> conn;
1318 for (auto sc : make_range(this->n_sub_elem()))
1319 {
1320 this->connectivity(sc, TECPLOT, conn);
1321
1322 std::copy(conn.begin(),
1323 conn.end(),
1324 std::ostream_iterator<dof_id_type>(out_stream, " "));
1325
1326 out_stream << '\n';
1327 }
1328 return;
1329 }
1330
1331 case UCD:
1332 {
1333 for (auto i : this->node_index_range())
1334 out_stream << this->node_id(i)+1 << "\t";
1335
1336 out_stream << '\n';
1337 return;
1338 }
1339
1340 default:
1341 libmesh_error_msg("Unsupported IO package " << iop);
1342 }
1343 }
1344
1345
1346
quality(const ElemQuality q)1347 Real Elem::quality (const ElemQuality q) const
1348 {
1349 switch (q)
1350 {
1351 /**
1352 * I don't know what to do for this metric.
1353 */
1354 default:
1355 {
1356 libmesh_do_once( libmesh_here();
1357
1358 libMesh::err << "ERROR: unknown quality metric: "
1359 << Utility::enum_to_string(q)
1360 << std::endl
1361 << "Cowardly returning 1."
1362 << std::endl; );
1363
1364 return 1.;
1365 }
1366 }
1367 }
1368
1369
1370
ancestor()1371 bool Elem::ancestor() const
1372 {
1373 #ifdef LIBMESH_ENABLE_AMR
1374
1375 // Use a fast, DistributedMesh-safe definition
1376 const bool is_ancestor =
1377 !this->active() && !this->subactive();
1378
1379 // But check for inconsistencies if we have time
1380 #ifdef DEBUG
1381 if (!is_ancestor && this->has_children())
1382 {
1383 for (auto & c : this->child_ref_range())
1384 {
1385 if (&c != remote_elem)
1386 {
1387 libmesh_assert(!c.active());
1388 libmesh_assert(!c.ancestor());
1389 }
1390 }
1391 }
1392 #endif // DEBUG
1393
1394 return is_ancestor;
1395
1396 #else
1397 return false;
1398 #endif
1399 }
1400
1401
1402
1403 #ifdef LIBMESH_ENABLE_AMR
1404
add_child(Elem * elem)1405 void Elem::add_child (Elem * elem)
1406 {
1407 const unsigned int nc = this->n_children();
1408
1409 if (_children == nullptr)
1410 {
1411 _children = new Elem *[nc];
1412
1413 for (unsigned int c = 0; c != nc; c++)
1414 this->set_child(c, nullptr);
1415 }
1416
1417 for (unsigned int c = 0; c != nc; c++)
1418 {
1419 if (this->_children[c] == nullptr || this->_children[c] == remote_elem)
1420 {
1421 libmesh_assert_equal_to (this, elem->parent());
1422 this->set_child(c, elem);
1423 return;
1424 }
1425 }
1426
1427 libmesh_error_msg("Error: Tried to add a child to an element with full children array");
1428 }
1429
1430
1431
add_child(Elem * elem,unsigned int c)1432 void Elem::add_child (Elem * elem, unsigned int c)
1433 {
1434 if (!this->has_children())
1435 {
1436 const unsigned int nc = this->n_children();
1437 _children = new Elem *[nc];
1438
1439 for (unsigned int i = 0; i != nc; i++)
1440 this->set_child(i, nullptr);
1441 }
1442
1443 libmesh_assert (this->_children[c] == nullptr || this->child_ptr(c) == remote_elem);
1444 libmesh_assert (elem == remote_elem || this == elem->parent());
1445
1446 this->set_child(c, elem);
1447 }
1448
1449
1450
replace_child(Elem * elem,unsigned int c)1451 void Elem::replace_child (Elem * elem, unsigned int c)
1452 {
1453 libmesh_assert(this->has_children());
1454
1455 libmesh_assert(this->child_ptr(c));
1456
1457 this->set_child(c, elem);
1458 }
1459
1460
1461
family_tree(std::vector<const Elem * > & family,bool reset)1462 void Elem::family_tree (std::vector<const Elem *> & family,
1463 bool reset) const
1464 {
1465 ElemInternal::family_tree(this, family, reset);
1466 }
1467
1468
1469
family_tree(std::vector<Elem * > & family,bool reset)1470 void Elem::family_tree (std::vector<Elem *> & family,
1471 bool reset)
1472 {
1473 ElemInternal::family_tree(this, family, reset);
1474 }
1475
1476
1477
total_family_tree(std::vector<const Elem * > & family,bool reset)1478 void Elem::total_family_tree (std::vector<const Elem *> & family,
1479 bool reset) const
1480 {
1481 ElemInternal::total_family_tree(this, family, reset);
1482 }
1483
1484
1485
total_family_tree(std::vector<Elem * > & family,bool reset)1486 void Elem::total_family_tree (std::vector<Elem *> & family,
1487 bool reset)
1488 {
1489 ElemInternal::total_family_tree(this, family, reset);
1490 }
1491
1492
1493
active_family_tree(std::vector<const Elem * > & active_family,bool reset)1494 void Elem::active_family_tree (std::vector<const Elem *> & active_family,
1495 bool reset) const
1496 {
1497 ElemInternal::active_family_tree(this, active_family, reset);
1498 }
1499
1500
1501
active_family_tree(std::vector<Elem * > & active_family,bool reset)1502 void Elem::active_family_tree (std::vector<Elem *> & active_family,
1503 bool reset)
1504 {
1505 ElemInternal::active_family_tree(this, active_family, reset);
1506 }
1507
1508
1509
family_tree_by_side(std::vector<const Elem * > & family,unsigned int side,bool reset)1510 void Elem::family_tree_by_side (std::vector<const Elem *> & family,
1511 unsigned int side,
1512 bool reset) const
1513 {
1514 ElemInternal::family_tree_by_side(this, family, side, reset);
1515 }
1516
1517
1518
family_tree_by_side(std::vector<Elem * > & family,unsigned int side,bool reset)1519 void Elem:: family_tree_by_side (std::vector<Elem *> & family,
1520 unsigned int side,
1521 bool reset)
1522 {
1523 ElemInternal::family_tree_by_side(this, family, side, reset);
1524 }
1525
1526
1527
active_family_tree_by_side(std::vector<const Elem * > & family,unsigned int side,bool reset)1528 void Elem::active_family_tree_by_side (std::vector<const Elem *> & family,
1529 unsigned int side,
1530 bool reset) const
1531 {
1532 ElemInternal::active_family_tree_by_side(this, family, side, reset);
1533 }
1534
1535
1536
active_family_tree_by_side(std::vector<Elem * > & family,unsigned int side,bool reset)1537 void Elem::active_family_tree_by_side (std::vector<Elem *> & family,
1538 unsigned int side,
1539 bool reset)
1540 {
1541 ElemInternal::active_family_tree_by_side(this, family, side, reset);
1542 }
1543
1544
1545
family_tree_by_neighbor(std::vector<const Elem * > & family,const Elem * neighbor,bool reset)1546 void Elem::family_tree_by_neighbor (std::vector<const Elem *> & family,
1547 const Elem * neighbor,
1548 bool reset) const
1549 {
1550 ElemInternal::family_tree_by_neighbor(this, family, neighbor, reset);
1551 }
1552
1553
1554
family_tree_by_neighbor(std::vector<Elem * > & family,Elem * neighbor,bool reset)1555 void Elem::family_tree_by_neighbor (std::vector<Elem *> & family,
1556 Elem * neighbor,
1557 bool reset)
1558 {
1559 ElemInternal::family_tree_by_neighbor(this, family, neighbor, reset);
1560 }
1561
1562
1563
total_family_tree_by_neighbor(std::vector<const Elem * > & family,const Elem * neighbor,bool reset)1564 void Elem::total_family_tree_by_neighbor (std::vector<const Elem *> & family,
1565 const Elem * neighbor,
1566 bool reset) const
1567 {
1568 ElemInternal::total_family_tree_by_neighbor(this, family, neighbor, reset);
1569 }
1570
1571
1572
total_family_tree_by_neighbor(std::vector<Elem * > & family,Elem * neighbor,bool reset)1573 void Elem::total_family_tree_by_neighbor (std::vector<Elem *> & family,
1574 Elem * neighbor,
1575 bool reset)
1576 {
1577 ElemInternal::total_family_tree_by_neighbor(this, family, neighbor, reset);
1578 }
1579
1580
1581
family_tree_by_subneighbor(std::vector<const Elem * > & family,const Elem * neighbor,const Elem * subneighbor,bool reset)1582 void Elem::family_tree_by_subneighbor (std::vector<const Elem *> & family,
1583 const Elem * neighbor,
1584 const Elem * subneighbor,
1585 bool reset) const
1586 {
1587 ElemInternal::family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
1588 }
1589
1590
1591
family_tree_by_subneighbor(std::vector<Elem * > & family,Elem * neighbor,Elem * subneighbor,bool reset)1592 void Elem::family_tree_by_subneighbor (std::vector<Elem *> & family,
1593 Elem * neighbor,
1594 Elem * subneighbor,
1595 bool reset)
1596 {
1597 ElemInternal::family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
1598 }
1599
1600
1601
total_family_tree_by_subneighbor(std::vector<const Elem * > & family,const Elem * neighbor,const Elem * subneighbor,bool reset)1602 void Elem::total_family_tree_by_subneighbor (std::vector<const Elem *> & family,
1603 const Elem * neighbor,
1604 const Elem * subneighbor,
1605 bool reset) const
1606 {
1607 ElemInternal::total_family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
1608 }
1609
1610
1611
total_family_tree_by_subneighbor(std::vector<Elem * > & family,Elem * neighbor,Elem * subneighbor,bool reset)1612 void Elem::total_family_tree_by_subneighbor (std::vector<Elem *> & family,
1613 Elem * neighbor,
1614 Elem * subneighbor,
1615 bool reset)
1616 {
1617 ElemInternal::total_family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
1618 }
1619
1620
1621
active_family_tree_by_neighbor(std::vector<const Elem * > & family,const Elem * neighbor,bool reset)1622 void Elem::active_family_tree_by_neighbor (std::vector<const Elem *> & family,
1623 const Elem * neighbor,
1624 bool reset) const
1625 {
1626 ElemInternal::active_family_tree_by_neighbor(this, family, neighbor, reset);
1627 }
1628
1629
1630
active_family_tree_by_neighbor(std::vector<Elem * > & family,Elem * neighbor,bool reset)1631 void Elem::active_family_tree_by_neighbor (std::vector<Elem *> & family,
1632 Elem * neighbor,
1633 bool reset)
1634 {
1635 ElemInternal::active_family_tree_by_neighbor(this, family, neighbor, reset);
1636 }
1637
1638
1639
active_family_tree_by_topological_neighbor(std::vector<const Elem * > & family,const Elem * neighbor,const MeshBase & mesh,const PointLocatorBase & point_locator,const PeriodicBoundaries * pb,bool reset)1640 void Elem::active_family_tree_by_topological_neighbor (std::vector<const Elem *> & family,
1641 const Elem * neighbor,
1642 const MeshBase & mesh,
1643 const PointLocatorBase & point_locator,
1644 const PeriodicBoundaries * pb,
1645 bool reset) const
1646 {
1647 ElemInternal::active_family_tree_by_topological_neighbor(this, family, neighbor,
1648 mesh, point_locator, pb,
1649 reset);
1650 }
1651
1652
1653
active_family_tree_by_topological_neighbor(std::vector<Elem * > & family,Elem * neighbor,const MeshBase & mesh,const PointLocatorBase & point_locator,const PeriodicBoundaries * pb,bool reset)1654 void Elem::active_family_tree_by_topological_neighbor (std::vector<Elem *> & family,
1655 Elem * neighbor,
1656 const MeshBase & mesh,
1657 const PointLocatorBase & point_locator,
1658 const PeriodicBoundaries * pb,
1659 bool reset)
1660 {
1661 ElemInternal::active_family_tree_by_topological_neighbor(this, family, neighbor,
1662 mesh, point_locator, pb,
1663 reset);
1664 }
1665
1666
is_child_on_edge(const unsigned int c,const unsigned int e)1667 bool Elem::is_child_on_edge(const unsigned int c,
1668 const unsigned int e) const
1669 {
1670 libmesh_assert_less (c, this->n_children());
1671 libmesh_assert_less (e, this->n_edges());
1672
1673 std::unique_ptr<const Elem> my_edge = this->build_edge_ptr(e);
1674 std::unique_ptr<const Elem> child_edge = this->child_ptr(c)->build_edge_ptr(e);
1675
1676 // We're assuming that an overlapping child edge has the same
1677 // number and orientation as its parent
1678 return (child_edge->node_id(0) == my_edge->node_id(0) ||
1679 child_edge->node_id(1) == my_edge->node_id(1));
1680 }
1681
1682
1683
min_p_level_by_neighbor(const Elem * neighbor_in,unsigned int current_min)1684 unsigned int Elem::min_p_level_by_neighbor(const Elem * neighbor_in,
1685 unsigned int current_min) const
1686 {
1687 libmesh_assert(!this->subactive());
1688 libmesh_assert(neighbor_in->active());
1689
1690 // If we're an active element this is simple
1691 if (this->active())
1692 return std::min(current_min, this->p_level());
1693
1694 libmesh_assert(has_neighbor(neighbor_in));
1695
1696 // The p_level() of an ancestor element is already the minimum
1697 // p_level() of its children - so if that's high enough, we don't
1698 // need to examine any children.
1699 if (current_min <= this->p_level())
1700 return current_min;
1701
1702 unsigned int min_p_level = current_min;
1703
1704 for (auto & c : this->child_ref_range())
1705 if (&c != remote_elem && c.has_neighbor(neighbor_in))
1706 min_p_level =
1707 c.min_p_level_by_neighbor(neighbor_in, min_p_level);
1708
1709 return min_p_level;
1710 }
1711
1712
min_new_p_level_by_neighbor(const Elem * neighbor_in,unsigned int current_min)1713 unsigned int Elem::min_new_p_level_by_neighbor(const Elem * neighbor_in,
1714 unsigned int current_min) const
1715 {
1716 libmesh_assert(!this->subactive());
1717 libmesh_assert(neighbor_in->active());
1718
1719 // If we're an active element this is simple
1720 if (this->active())
1721 {
1722 unsigned int new_p_level = this->p_level();
1723 if (this->p_refinement_flag() == Elem::REFINE)
1724 new_p_level += 1;
1725 if (this->p_refinement_flag() == Elem::COARSEN)
1726 {
1727 libmesh_assert_greater (new_p_level, 0);
1728 new_p_level -= 1;
1729 }
1730 return std::min(current_min, new_p_level);
1731 }
1732
1733 libmesh_assert(has_neighbor(neighbor_in));
1734
1735 unsigned int min_p_level = current_min;
1736
1737 for (auto & c : this->child_ref_range())
1738 if (&c != remote_elem && c.has_neighbor(neighbor_in))
1739 min_p_level =
1740 c.min_new_p_level_by_neighbor(neighbor_in, min_p_level);
1741
1742 return min_p_level;
1743 }
1744
1745
1746
as_parent_node(unsigned int child,unsigned int child_node)1747 unsigned int Elem::as_parent_node (unsigned int child,
1748 unsigned int child_node) const
1749 {
1750 const unsigned int nc = this->n_children();
1751 libmesh_assert_less(child, nc);
1752
1753 // Cached return values, indexed first by embedding_matrix version,
1754 // then by child number, then by child node number.
1755 std::vector<std::vector<std::vector<signed char>>> &
1756 cached_parent_indices = this->_get_parent_indices_cache();
1757
1758 unsigned int em_vers = this->embedding_matrix_version();
1759
1760 // We may be updating the cache on one thread, and while that
1761 // happens we can't safely access the cache from other threads.
1762 Threads::spin_mutex::scoped_lock lock(parent_indices_mutex);
1763
1764 if (em_vers >= cached_parent_indices.size())
1765 cached_parent_indices.resize(em_vers+1);
1766
1767 if (child >= cached_parent_indices[em_vers].size())
1768 {
1769 const signed char nn = cast_int<signed char>(this->n_nodes());
1770
1771 cached_parent_indices[em_vers].resize(nc);
1772
1773 for (unsigned int c = 0; c != nc; ++c)
1774 {
1775 const unsigned int ncn = this->n_nodes_in_child(c);
1776 cached_parent_indices[em_vers][c].resize(ncn);
1777 for (unsigned int cn = 0; cn != ncn; ++cn)
1778 {
1779 for (signed char n = 0; n != nn; ++n)
1780 {
1781 const float em_val = this->embedding_matrix
1782 (c, cn, n);
1783 if (em_val == 1)
1784 {
1785 cached_parent_indices[em_vers][c][cn] = n;
1786 break;
1787 }
1788
1789 if (em_val != 0)
1790 {
1791 cached_parent_indices[em_vers][c][cn] =
1792 -1;
1793 break;
1794 }
1795
1796 // We should never see an all-zero embedding matrix
1797 // row
1798 libmesh_assert_not_equal_to (n+1, nn);
1799 }
1800 }
1801 }
1802 }
1803
1804 const signed char cache_val =
1805 cached_parent_indices[em_vers][child][child_node];
1806 if (cache_val == -1)
1807 return libMesh::invalid_uint;
1808
1809 return cached_parent_indices[em_vers][child][child_node];
1810 }
1811
1812
1813
1814 const std::vector<std::pair<unsigned char, unsigned char>> &
parent_bracketing_nodes(unsigned int child,unsigned int child_node)1815 Elem::parent_bracketing_nodes(unsigned int child,
1816 unsigned int child_node) const
1817 {
1818 // Indexed first by embedding matrix type, then by child id, then by
1819 // child node, then by bracketing pair
1820 std::vector<std::vector<std::vector<std::vector<std::pair<unsigned char, unsigned char>>>>> &
1821 cached_bracketing_nodes = this->_get_bracketing_node_cache();
1822
1823 const unsigned int em_vers = this->embedding_matrix_version();
1824
1825 // We may be updating the cache on one thread, and while that
1826 // happens we can't safely access the cache from other threads.
1827 Threads::spin_mutex::scoped_lock lock(parent_bracketing_nodes_mutex);
1828
1829 if (cached_bracketing_nodes.size() <= em_vers)
1830 cached_bracketing_nodes.resize(em_vers+1);
1831
1832 const unsigned int nc = this->n_children();
1833
1834 // If we haven't cached the bracketing nodes corresponding to this
1835 // embedding matrix yet, let's do so now.
1836 if (cached_bracketing_nodes[em_vers].size() < nc)
1837 {
1838 // If we're a second-order element but we're not a full-order
1839 // element, then some of our bracketing nodes may not exist
1840 // except on the equivalent full-order element. Let's build an
1841 // equivalent full-order element and make a copy of its cache to
1842 // use.
1843 if (this->default_order() != FIRST &&
1844 second_order_equivalent_type(this->type(), /*full_ordered=*/ true) != this->type())
1845 {
1846 // Check that we really are the non-full-order type
1847 libmesh_assert_equal_to
1848 (second_order_equivalent_type (this->type(), false),
1849 this->type());
1850
1851 // Build the full-order type
1852 ElemType full_type =
1853 second_order_equivalent_type(this->type(), /*full_ordered=*/ true);
1854 std::unique_ptr<Elem> full_elem = Elem::build(full_type);
1855
1856 // This won't work for elements with multiple
1857 // embedding_matrix versions, but every such element is full
1858 // order anyways.
1859 libmesh_assert_equal_to(em_vers, 0);
1860
1861 // Make sure its cache has been built. We temporarily
1862 // release our mutex lock so that the inner call can
1863 // re-acquire it.
1864 lock.release();
1865 full_elem->parent_bracketing_nodes(0,0);
1866
1867 // And then we need to lock again, so that if someone *else*
1868 // grabbed our lock before we did we don't risk accessing
1869 // cached_bracketing_nodes while they're working on it.
1870 // Threading is hard.
1871 lock.acquire(parent_bracketing_nodes_mutex);
1872
1873 // Copy its cache
1874 cached_bracketing_nodes =
1875 full_elem->_get_bracketing_node_cache();
1876
1877 // Now we don't need to build the cache ourselves.
1878 return cached_bracketing_nodes[em_vers][child][child_node];
1879 }
1880
1881 cached_bracketing_nodes[em_vers].resize(nc);
1882
1883 const unsigned int nn = this->n_nodes();
1884
1885 // We have to examine each child
1886 for (unsigned int c = 0; c != nc; ++c)
1887 {
1888 const unsigned int ncn = this->n_nodes_in_child(c);
1889
1890 cached_bracketing_nodes[em_vers][c].resize(ncn);
1891
1892 // We have to examine each node in that child
1893 for (unsigned int n = 0; n != ncn; ++n)
1894 {
1895 // If this child node isn't a vertex or an infinite
1896 // child element's mid-infinite-edge node, then we need
1897 // to find bracketing nodes on the child.
1898 if (!this->is_vertex_on_child(c, n)
1899 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1900 && !this->is_mid_infinite_edge_node(n)
1901 #endif
1902 )
1903 {
1904 // Use the embedding matrix to find the child node
1905 // location in parent master element space
1906 Point bracketed_pt;
1907
1908 for (unsigned int pn = 0; pn != nn; ++pn)
1909 {
1910 const float em_val =
1911 this->embedding_matrix(c,n,pn);
1912
1913 libmesh_assert_not_equal_to (em_val, 1);
1914 if (em_val != 0.)
1915 bracketed_pt.add_scaled(this->master_point(pn), em_val);
1916 }
1917
1918 // Check each pair of nodes on the child which are
1919 // also both parent nodes
1920 for (unsigned int n1 = 0; n1 != ncn; ++n1)
1921 {
1922 if (n1 == n)
1923 continue;
1924
1925 unsigned int parent_n1 =
1926 this->as_parent_node(c,n1);
1927
1928 if (parent_n1 == libMesh::invalid_uint)
1929 continue;
1930
1931 Point p1 = this->master_point(parent_n1);
1932
1933 for (unsigned int n2 = n1+1; n2 < nn; ++n2)
1934 {
1935 if (n2 == n)
1936 continue;
1937
1938 unsigned int parent_n2 =
1939 this->as_parent_node(c,n2);
1940
1941 if (parent_n2 == libMesh::invalid_uint)
1942 continue;
1943
1944 Point p2 = this->master_point(parent_n2);
1945
1946 Point pmid = (p1 + p2)/2;
1947
1948 if (pmid == bracketed_pt)
1949 {
1950 cached_bracketing_nodes[em_vers][c][n].emplace_back(parent_n1, parent_n2);
1951 break;
1952 }
1953 else
1954 libmesh_assert(!pmid.absolute_fuzzy_equals(bracketed_pt));
1955 }
1956 }
1957 }
1958 // If this child node is a parent node, we need to
1959 // find bracketing nodes on the parent.
1960 else
1961 {
1962 unsigned int parent_node = this->as_parent_node(c,n);
1963
1964 Point bracketed_pt;
1965
1966 // If we're not a parent node, use the embedding
1967 // matrix to find the child node location in parent
1968 // master element space
1969 if (parent_node == libMesh::invalid_uint)
1970 {
1971 for (unsigned int pn = 0; pn != nn; ++pn)
1972 {
1973 const float em_val =
1974 this->embedding_matrix(c,n,pn);
1975
1976 libmesh_assert_not_equal_to (em_val, 1);
1977 if (em_val != 0.)
1978 bracketed_pt.add_scaled(this->master_point(pn), em_val);
1979 }
1980 }
1981 // If we're a parent node then we need no arithmetic
1982 else
1983 bracketed_pt = this->master_point(parent_node);
1984
1985 for (unsigned int n1 = 0; n1 != nn; ++n1)
1986 {
1987 if (n1 == parent_node)
1988 continue;
1989
1990 Point p1 = this->master_point(n1);
1991
1992 for (unsigned int n2 = n1+1; n2 < nn; ++n2)
1993 {
1994 if (n2 == parent_node)
1995 continue;
1996
1997 Point pmid = (p1 + this->master_point(n2))/2;
1998
1999 if (pmid == bracketed_pt)
2000 {
2001 cached_bracketing_nodes[em_vers][c][n].emplace_back(n1, n2);
2002 break;
2003 }
2004 else
2005 libmesh_assert(!pmid.absolute_fuzzy_equals(bracketed_pt));
2006 }
2007 }
2008 }
2009 }
2010 }
2011 }
2012
2013 return cached_bracketing_nodes[em_vers][child][child_node];
2014 }
2015
2016
2017 const std::vector<std::pair<dof_id_type, dof_id_type>>
bracketing_nodes(unsigned int child,unsigned int child_node)2018 Elem::bracketing_nodes(unsigned int child,
2019 unsigned int child_node) const
2020 {
2021 std::vector<std::pair<dof_id_type, dof_id_type>> returnval;
2022
2023 const std::vector<std::pair<unsigned char, unsigned char>> & pbc =
2024 this->parent_bracketing_nodes(child,child_node);
2025
2026 for (const auto & pb : pbc)
2027 {
2028 const unsigned short n_n = this->n_nodes();
2029 if (pb.first < n_n && pb.second < n_n)
2030 returnval.emplace_back(this->node_id(pb.first), this->node_id(pb.second));
2031 else
2032 {
2033 // We must be on a non-full-order higher order element...
2034 libmesh_assert_not_equal_to(this->default_order(), FIRST);
2035 libmesh_assert_not_equal_to
2036 (second_order_equivalent_type (this->type(), true),
2037 this->type());
2038 libmesh_assert_equal_to
2039 (second_order_equivalent_type (this->type(), false),
2040 this->type());
2041
2042 // And that's a shame, because this is a nasty search:
2043
2044 // Build the full-order type
2045 ElemType full_type =
2046 second_order_equivalent_type(this->type(), /*full_ordered=*/ true);
2047 std::unique_ptr<Elem> full_elem = Elem::build(full_type);
2048
2049 dof_id_type pt1 = DofObject::invalid_id;
2050 dof_id_type pt2 = DofObject::invalid_id;
2051
2052 // Find the bracketing nodes by figuring out what
2053 // already-created children will have them.
2054
2055 // This only doesn't break horribly because we add children
2056 // and nodes in straightforward + hierarchical orders...
2057 for (unsigned int c=0; c <= child; ++c)
2058 for (auto n : make_range(this->n_nodes_in_child(c)))
2059 {
2060 if (c == child && n == child_node)
2061 break;
2062
2063 if (pb.first == full_elem->as_parent_node(c,n))
2064 {
2065 // We should be consistent
2066 if (pt1 != DofObject::invalid_id)
2067 libmesh_assert_equal_to(pt1, this->child_ptr(c)->node_id(n));
2068
2069 pt1 = this->child_ptr(c)->node_id(n);
2070 }
2071
2072 if (pb.second == full_elem->as_parent_node(c,n))
2073 {
2074 // We should be consistent
2075 if (pt2 != DofObject::invalid_id)
2076 libmesh_assert_equal_to(pt2, this->child_ptr(c)->node_id(n));
2077
2078 pt2 = this->child_ptr(c)->node_id(n);
2079 }
2080 }
2081
2082 // We should *usually* find all bracketing nodes by the time
2083 // we query them (again, because of the child & node add
2084 // order)
2085 //
2086 // The exception is if we're a HEX20, in which case we will
2087 // find pairs of vertex nodes and edge nodes bracketing the
2088 // new central node but we *won't* find the pairs of face
2089 // nodes which we would have had on a HEX27. In that case
2090 // we'll still have enough bracketing nodes for a
2091 // topological lookup, but we won't be able to make the
2092 // following assertions.
2093 if (this->type() != HEX20)
2094 {
2095 libmesh_assert_not_equal_to (pt1, DofObject::invalid_id);
2096 libmesh_assert_not_equal_to (pt2, DofObject::invalid_id);
2097 }
2098
2099 if (pt1 != DofObject::invalid_id &&
2100 pt2 != DofObject::invalid_id)
2101 returnval.emplace_back(pt1, pt2);
2102 }
2103 }
2104
2105 return returnval;
2106 }
2107 #endif // #ifdef LIBMESH_ENABLE_AMR
2108
2109
2110
2111
contains_point(const Point & p,Real tol)2112 bool Elem::contains_point (const Point & p, Real tol) const
2113 {
2114 // We currently allow the user to enlarge the bounding box by
2115 // providing a tol > TOLERANCE (so this routine is identical to
2116 // Elem::close_to_point()), but print a warning so that the
2117 // user can eventually switch his code over to calling close_to_point()
2118 // instead, which is intended to be used for this purpose.
2119 if (tol > TOLERANCE)
2120 {
2121 libmesh_do_once(libMesh::err
2122 << "WARNING: Resizing bounding box to match user-specified tolerance!\n"
2123 << "In the future, calls to Elem::contains_point() with tol > TOLERANCE\n"
2124 << "will be more optimized, but should not be used\n"
2125 << "to search for points 'close to' elements!\n"
2126 << "Instead, use Elem::close_to_point() for this purpose.\n"
2127 << std::endl;);
2128 return this->point_test(p, tol, tol);
2129 }
2130 else
2131 return this->point_test(p, TOLERANCE, tol);
2132 }
2133
2134
2135
2136
close_to_point(const Point & p,Real tol)2137 bool Elem::close_to_point (const Point & p, Real tol) const
2138 {
2139 // This test uses the user's passed-in tolerance for the
2140 // bounding box test as well, thereby allowing the routine to
2141 // find points which are not only "in" the element, but also
2142 // "nearby" to within some tolerance.
2143 return this->point_test(p, tol, tol);
2144 }
2145
2146
2147
2148
point_test(const Point & p,Real box_tol,Real map_tol)2149 bool Elem::point_test(const Point & p, Real box_tol, Real map_tol) const
2150 {
2151 libmesh_assert_greater (box_tol, 0.);
2152 libmesh_assert_greater (map_tol, 0.);
2153
2154 // This is a great optimization on first order elements, but it
2155 // could return false negatives on higher orders
2156 if (this->default_order() == FIRST)
2157 {
2158 // Check to make sure the element *could* contain this point, so we
2159 // can avoid an expensive inverse_map call if it doesn't.
2160 bool
2161 #if LIBMESH_DIM > 2
2162 point_above_min_z = false,
2163 point_below_max_z = false,
2164 #endif
2165 #if LIBMESH_DIM > 1
2166 point_above_min_y = false,
2167 point_below_max_y = false,
2168 #endif
2169 point_above_min_x = false,
2170 point_below_max_x = false;
2171
2172 // For relative bounding box checks in physical space
2173 const Real my_hmax = this->hmax();
2174
2175 for (auto & n : this->node_ref_range())
2176 {
2177 point_above_min_x = point_above_min_x || (n(0) - my_hmax*box_tol <= p(0));
2178 point_below_max_x = point_below_max_x || (n(0) + my_hmax*box_tol >= p(0));
2179 #if LIBMESH_DIM > 1
2180 point_above_min_y = point_above_min_y || (n(1) - my_hmax*box_tol <= p(1));
2181 point_below_max_y = point_below_max_y || (n(1) + my_hmax*box_tol >= p(1));
2182 #endif
2183 #if LIBMESH_DIM > 2
2184 point_above_min_z = point_above_min_z || (n(2) - my_hmax*box_tol <= p(2));
2185 point_below_max_z = point_below_max_z || (n(2) + my_hmax*box_tol >= p(2));
2186 #endif
2187 }
2188
2189 if (
2190 #if LIBMESH_DIM > 2
2191 !point_above_min_z ||
2192 !point_below_max_z ||
2193 #endif
2194 #if LIBMESH_DIM > 1
2195 !point_above_min_y ||
2196 !point_below_max_y ||
2197 #endif
2198 !point_above_min_x ||
2199 !point_below_max_x)
2200 return false;
2201 }
2202
2203 // To be on the safe side, we converge the inverse_map() iteration
2204 // to a slightly tighter tolerance than that requested by the
2205 // user...
2206 const Point mapped_point = FEMap::inverse_map(this->dim(),
2207 this,
2208 p,
2209 0.1*map_tol, // <- this is |dx| tolerance, the Newton residual should be ~ |dx|^2
2210 /*secure=*/ false,
2211 /*extra_checks=*/ false);
2212
2213 // Check that the refspace point maps back to p! This is only necessary
2214 // for 1D and 2D elements, 3D elements always live in 3D.
2215 //
2216 // TODO: The contains_point() function could most likely be implemented
2217 // more efficiently in the element sub-classes themselves, at least for
2218 // the linear element types.
2219 if (this->dim() < 3)
2220 {
2221 Point xyz = FEMap::map(this->dim(), this, mapped_point);
2222
2223 // Compute the distance between the original point and the re-mapped point.
2224 // They should be in the same place.
2225 Real dist = (xyz - p).norm();
2226
2227
2228 // If dist is larger than some fraction of the tolerance, then return false.
2229 // This can happen when e.g. a 2D element is living in 3D, and
2230 // FEMap::inverse_map() maps p onto the projection of the element,
2231 // effectively "tricking" FEInterface::on_reference_element().
2232 if (dist > this->hmax() * map_tol)
2233 return false;
2234 }
2235
2236
2237
2238 return FEInterface::on_reference_element(mapped_point, this->type(), map_tol);
2239 }
2240
2241
2242
2243
print_info(std::ostream & os)2244 void Elem::print_info (std::ostream & os) const
2245 {
2246 os << this->get_info()
2247 << std::endl;
2248 }
2249
2250
2251
get_info()2252 std::string Elem::get_info () const
2253 {
2254 std::ostringstream oss;
2255
2256 oss << " Elem Information" << '\n'
2257 << " id()=";
2258
2259 if (this->valid_id())
2260 oss << this->id();
2261 else
2262 oss << "invalid";
2263
2264 #ifdef LIBMESH_ENABLE_UNIQUE_ID
2265 oss << ", unique_id()=";
2266 if (this->valid_unique_id())
2267 oss << this->unique_id();
2268 else
2269 oss << "invalid";
2270 #endif
2271
2272 oss << ", processor_id()=" << this->processor_id() << '\n';
2273
2274 oss << " type()=" << Utility::enum_to_string(this->type()) << '\n'
2275 << " dim()=" << this->dim() << '\n'
2276 << " n_nodes()=" << this->n_nodes() << '\n';
2277
2278 for (auto n : this->node_index_range())
2279 oss << " " << n << this->node_ref(n);
2280
2281 oss << " n_sides()=" << this->n_sides() << '\n';
2282
2283 for (auto s : this->side_index_range())
2284 {
2285 oss << " neighbor(" << s << ")=";
2286 if (this->neighbor_ptr(s))
2287 oss << this->neighbor_ptr(s)->id() << '\n';
2288 else
2289 oss << "nullptr\n";
2290 }
2291
2292 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2293 if (!this->infinite())
2294 {
2295 #endif
2296 oss << " hmin()=" << this->hmin()
2297 << ", hmax()=" << this->hmax() << '\n'
2298 << " volume()=" << this->volume() << '\n';
2299 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2300 }
2301 #endif
2302 oss << " active()=" << this->active()
2303 << ", ancestor()=" << this->ancestor()
2304 << ", subactive()=" << this->subactive()
2305 << ", has_children()=" << this->has_children() << '\n'
2306 << " parent()=";
2307 if (this->parent())
2308 oss << this->parent()->id() << '\n';
2309 else
2310 oss << "nullptr\n";
2311 oss << " level()=" << this->level()
2312 << ", p_level()=" << this->p_level() << '\n'
2313 #ifdef LIBMESH_ENABLE_AMR
2314 << " refinement_flag()=" << Utility::enum_to_string(this->refinement_flag()) << '\n'
2315 << " p_refinement_flag()=" << Utility::enum_to_string(this->p_refinement_flag()) << '\n'
2316 #endif
2317 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2318 << " infinite()=" << this->infinite() << '\n';
2319 if (this->infinite())
2320 oss << " origin()=" << this->origin() << '\n'
2321 #endif
2322 ;
2323
2324 oss << " DoFs=";
2325 for (auto s : make_range(this->n_systems()))
2326 for (auto v : make_range(this->n_vars(s)))
2327 for (auto c : make_range(this->n_comp(s,v)))
2328 oss << '(' << s << '/' << v << '/' << this->dof_number(s,v,c) << ") ";
2329
2330
2331 return oss.str();
2332 }
2333
2334
2335
nullify_neighbors()2336 void Elem::nullify_neighbors ()
2337 {
2338 // Tell any of my neighbors about my death...
2339 // Looks strange, huh?
2340 for (auto n : this->side_index_range())
2341 {
2342 Elem * current_neighbor = this->neighbor_ptr(n);
2343 if (current_neighbor && current_neighbor != remote_elem)
2344 {
2345 // Note: it is possible that I see the neighbor
2346 // (which is coarser than me)
2347 // but they don't see me, so avoid that case.
2348 if (current_neighbor->level() == this->level())
2349 {
2350 const unsigned int w_n_a_i = current_neighbor->which_neighbor_am_i(this);
2351 libmesh_assert_less (w_n_a_i, current_neighbor->n_neighbors());
2352 current_neighbor->set_neighbor(w_n_a_i, nullptr);
2353 this->set_neighbor(n, nullptr);
2354 }
2355 }
2356 }
2357 }
2358
2359
2360
2361
n_second_order_adjacent_vertices(const unsigned int)2362 unsigned int Elem::n_second_order_adjacent_vertices (const unsigned int) const
2363 {
2364 // for linear elements, always return 0
2365 return 0;
2366 }
2367
2368
2369
second_order_adjacent_vertex(const unsigned int,const unsigned int)2370 unsigned short int Elem::second_order_adjacent_vertex (const unsigned int,
2371 const unsigned int) const
2372 {
2373 // for linear elements, always return 0
2374 return 0;
2375 }
2376
2377
2378
2379 std::pair<unsigned short int, unsigned short int>
second_order_child_vertex(const unsigned int)2380 Elem::second_order_child_vertex (const unsigned int) const
2381 {
2382 // for linear elements, always return 0
2383 return std::pair<unsigned short int, unsigned short int>(0,0);
2384 }
2385
2386
2387
first_order_equivalent_type(const ElemType et)2388 ElemType Elem::first_order_equivalent_type (const ElemType et)
2389 {
2390 switch (et)
2391 {
2392 case NODEELEM:
2393 return NODEELEM;
2394 case EDGE2:
2395 case EDGE3:
2396 case EDGE4:
2397 return EDGE2;
2398 case TRI3:
2399 case TRI6:
2400 return TRI3;
2401 case TRISHELL3:
2402 return TRISHELL3;
2403 case QUAD4:
2404 case QUAD8:
2405 case QUAD9:
2406 return QUAD4;
2407 case QUADSHELL4:
2408 case QUADSHELL8:
2409 return QUADSHELL4;
2410 case TET4:
2411 case TET10:
2412 return TET4;
2413 case HEX8:
2414 case HEX27:
2415 case HEX20:
2416 return HEX8;
2417 case PRISM6:
2418 case PRISM15:
2419 case PRISM18:
2420 return PRISM6;
2421 case PYRAMID5:
2422 case PYRAMID13:
2423 case PYRAMID14:
2424 return PYRAMID5;
2425
2426 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2427
2428 case INFEDGE2:
2429 return INFEDGE2;
2430 case INFQUAD4:
2431 case INFQUAD6:
2432 return INFQUAD4;
2433 case INFHEX8:
2434 case INFHEX16:
2435 case INFHEX18:
2436 return INFHEX8;
2437 case INFPRISM6:
2438 case INFPRISM12:
2439 return INFPRISM6;
2440
2441 #endif
2442
2443 default:
2444 // unknown element
2445 return INVALID_ELEM;
2446 }
2447 }
2448
2449
2450
second_order_equivalent_type(const ElemType et,const bool full_ordered)2451 ElemType Elem::second_order_equivalent_type (const ElemType et,
2452 const bool full_ordered)
2453 {
2454 switch (et)
2455 {
2456 case NODEELEM:
2457 return NODEELEM;
2458 case EDGE2:
2459 case EDGE3:
2460 {
2461 // full_ordered not relevant
2462 return EDGE3;
2463 }
2464
2465 case EDGE4:
2466 {
2467 // full_ordered not relevant
2468 return EDGE4;
2469 }
2470
2471 case TRI3:
2472 case TRI6:
2473 {
2474 // full_ordered not relevant
2475 return TRI6;
2476 }
2477
2478 // Currently there is no TRISHELL6, so similarly to other types
2479 // where this is the case, we just return the input.
2480 case TRISHELL3:
2481 return TRISHELL3;
2482
2483 case QUAD4:
2484 case QUAD8:
2485 {
2486 if (full_ordered)
2487 return QUAD9;
2488 else
2489 return QUAD8;
2490 }
2491
2492 case QUADSHELL4:
2493 case QUADSHELL8:
2494 {
2495 // There is no QUADSHELL9, so in that sense QUADSHELL8 is the
2496 // "full ordered" element.
2497 return QUADSHELL8;
2498 }
2499
2500 case QUAD9:
2501 {
2502 // full_ordered not relevant
2503 return QUAD9;
2504 }
2505
2506 case TET4:
2507 case TET10:
2508 {
2509 // full_ordered not relevant
2510 return TET10;
2511 }
2512
2513 case HEX8:
2514 case HEX20:
2515 {
2516 // see below how this correlates with INFHEX8
2517 if (full_ordered)
2518 return HEX27;
2519 else
2520 return HEX20;
2521 }
2522
2523 case HEX27:
2524 {
2525 // full_ordered not relevant
2526 return HEX27;
2527 }
2528
2529 case PRISM6:
2530 case PRISM15:
2531 {
2532 if (full_ordered)
2533 return PRISM18;
2534 else
2535 return PRISM15;
2536 }
2537
2538 case PRISM18:
2539 {
2540 // full_ordered not relevant
2541 return PRISM18;
2542 }
2543
2544 case PYRAMID5:
2545 case PYRAMID13:
2546 {
2547 if (full_ordered)
2548 return PYRAMID14;
2549 else
2550 return PYRAMID13;
2551 }
2552
2553 case PYRAMID14:
2554 {
2555 // full_ordered not relevant
2556 return PYRAMID14;
2557 }
2558
2559
2560
2561 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2562
2563 // infinite elements
2564 case INFEDGE2:
2565 {
2566 return INFEDGE2;
2567 }
2568
2569 case INFQUAD4:
2570 case INFQUAD6:
2571 {
2572 // full_ordered not relevant
2573 return INFQUAD6;
2574 }
2575
2576 case INFHEX8:
2577 case INFHEX16:
2578 {
2579 /*
2580 * Note that this matches with \p Hex8:
2581 * For full-ordered, \p InfHex18 and \p Hex27
2582 * belong together, and for not full-ordered,
2583 * \p InfHex16 and \p Hex20 belong together.
2584 */
2585 if (full_ordered)
2586 return INFHEX18;
2587 else
2588 return INFHEX16;
2589 }
2590
2591 case INFHEX18:
2592 {
2593 // full_ordered not relevant
2594 return INFHEX18;
2595 }
2596
2597 case INFPRISM6:
2598 case INFPRISM12:
2599 {
2600 // full_ordered not relevant
2601 return INFPRISM12;
2602 }
2603
2604 #endif
2605
2606
2607 default:
2608 {
2609 // what did we miss?
2610 libmesh_error_msg("No second order equivalent element type for et = "
2611 << Utility::enum_to_string(et));
2612 }
2613 }
2614 }
2615
2616
2617
boundary_sides_begin()2618 Elem::side_iterator Elem::boundary_sides_begin()
2619 {
2620 Predicates::BoundarySide<SideIter> bsp;
2621 return side_iterator(this->_first_side(), this->_last_side(), bsp);
2622 }
2623
2624
2625
2626
boundary_sides_end()2627 Elem::side_iterator Elem::boundary_sides_end()
2628 {
2629 Predicates::BoundarySide<SideIter> bsp;
2630 return side_iterator(this->_last_side(), this->_last_side(), bsp);
2631 }
2632
2633
2634
2635
volume()2636 Real Elem::volume () const
2637 {
2638 // The default implementation builds a finite element of the correct
2639 // order and sums up the JxW contributions. This can be expensive,
2640 // so the various element types can overload this method and compute
2641 // the volume more efficiently.
2642 const FEFamily mapping_family = FEMap::map_fe_type(*this);
2643 const FEType fe_type(this->default_order(), mapping_family);
2644
2645 std::unique_ptr<FEBase> fe (FEBase::build(this->dim(),
2646 fe_type));
2647
2648 const std::vector<Real> & JxW = fe->get_JxW();
2649
2650 // The default quadrature rule should integrate the mass matrix,
2651 // thus it should be plenty to compute the area
2652 QGauss qrule (this->dim(), fe_type.default_quadrature_order());
2653
2654 fe->attach_quadrature_rule(&qrule);
2655
2656 fe->reinit(this);
2657
2658 Real vol=0.;
2659 for (auto jxw : JxW)
2660 vol += jxw;
2661
2662 return vol;
2663
2664 }
2665
2666
2667
loose_bounding_box()2668 BoundingBox Elem::loose_bounding_box () const
2669 {
2670 Point pmin = this->point(0);
2671 Point pmax = pmin;
2672
2673 unsigned int n_points = this->n_nodes();
2674 for (unsigned int p=0; p != n_points; ++p)
2675 for (unsigned d=0; d<LIBMESH_DIM; ++d)
2676 {
2677 const Point & pt = this->point(p);
2678 if (pmin(d) > pt(d))
2679 pmin(d) = pt(d);
2680
2681 if (pmax(d) < pt(d))
2682 pmax(d) = pt(d);
2683 }
2684
2685 return BoundingBox(pmin, pmax);
2686 }
2687
2688
2689
is_vertex_on_parent(unsigned int c,unsigned int n)2690 bool Elem::is_vertex_on_parent(unsigned int c,
2691 unsigned int n) const
2692 {
2693 #ifdef LIBMESH_ENABLE_AMR
2694
2695 unsigned int my_n_vertices = this->n_vertices();
2696 for (unsigned int n_parent = 0; n_parent != my_n_vertices;
2697 ++n_parent)
2698 if (this->node_ptr(n_parent) == this->child_ptr(c)->node_ptr(n))
2699 return true;
2700 return false;
2701
2702 #else
2703
2704 // No AMR?
2705 libmesh_ignore(c,n);
2706 libmesh_error_msg("ERROR: AMR disabled, how did we get here?");
2707 return true;
2708
2709 #endif
2710 }
2711
2712
2713
opposite_side(const unsigned int)2714 unsigned int Elem::opposite_side(const unsigned int /*s*/) const
2715 {
2716 // If the subclass didn't rederive this, using it is an error
2717 libmesh_not_implemented();
2718 }
2719
2720
2721
opposite_node(const unsigned int,const unsigned int)2722 unsigned int Elem::opposite_node(const unsigned int /*n*/,
2723 const unsigned int /*s*/) const
2724 {
2725 // If the subclass didn't rederive this, using it is an error
2726 libmesh_not_implemented();
2727 }
2728
2729 } // namespace libMesh
2730