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