LCOV - code coverage report
Current view: top level - include/geom - elem.h (source / functions) Hit Total Coverage
Test: libMesh/libmesh: #4476 (4beb67) with base a68cc6 Lines: 355 436 81.4 %
Date: 2026-06-03 20:22:46 Functions: 142 182 78.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : 
      20             : #ifndef LIBMESH_ELEM_H
      21             : #define LIBMESH_ELEM_H
      22             : 
      23             : // Local includes
      24             : #include "libmesh/libmesh_common.h"
      25             : #include "libmesh/bounding_box.h"
      26             : #include "libmesh/dof_object.h"
      27             : #include "libmesh/id_types.h"
      28             : #include "libmesh/reference_counted_object.h"
      29             : #include "libmesh/node.h"
      30             : #include "libmesh/enum_elem_type.h" // INVALID_ELEM
      31             : #include "libmesh/multi_predicates.h"
      32             : #include "libmesh/pointer_to_pointer_iter.h"
      33             : #include "libmesh/int_range.h"
      34             : #include "libmesh/simple_range.h"
      35             : #include "libmesh/variant_filter_iterator.h"
      36             : #include "libmesh/hashword.h" // Used in compute_key() functions
      37             : 
      38             : // C++ includes
      39             : #include <algorithm>
      40             : #include <cstddef>
      41             : #include <iostream>
      42             : #include <limits.h> // CHAR_BIT
      43             : #include <set>
      44             : #include <vector>
      45             : #include <memory>
      46             : #include <array>
      47             : 
      48             : namespace libMesh
      49             : {
      50             : 
      51             : // Forward declarations
      52             : class BoundaryInfo;
      53             : class Elem;
      54             : class MeshBase;
      55             : class MeshRefinement;
      56             : #ifdef LIBMESH_ENABLE_PERIODIC
      57             : class PeriodicBoundaries;
      58             : class PointLocatorBase;
      59             : #endif
      60             : template <class SideType, class ParentType>
      61             : class Side;
      62             : enum ElemQuality : int;
      63             : enum IOPackage : int;
      64             : enum Order : int;
      65             : 
      66             : 
      67             : /**
      68             :  * This is the base class from which all geometric element types are
      69             :  * derived.  The \p Elem class provides standard information such as
      70             :  * the number of nodes, edges, faces, vertices, children, and
      71             :  * neighbors it has, as well as access to (or the ability to
      72             :  * construct) these entities.
      73             :  *
      74             :  * An \p Elem has pointers to its \p Node objects.  Some of these
      75             :  * nodes live at the vertices of the element, while others may live on
      76             :  * edges (and faces in 3D), or interior to the element.  The number of
      77             :  * nodes in a given element, \p n_nodes(), is encoded into the name of
      78             :  * the class.  For example, a \p Tri3 has three nodes which correspond
      79             :  * to the vertices, while a \p Tri6 has six nodes, three of which are
      80             :  * located at vertices, and three which are located at the midpoint of
      81             :  * each edge.  Nodes on edges, faces, and element interiors are called
      82             :  * second-order nodes.
      83             :  *
      84             :  * A 1D Elem is an \p Edge, a 2D Elem is a \p Face, and a 3D Elem is a
      85             :  * \p Cell.  An \p Elem is composed of a number of sides, which can
      86             :  * be accessed as dim-1 dimensional \p Elem types.  For example, a \p
      87             :  * Hex8 is a 3D hexahedral element. A \p Hex8 has 6 sides, which are
      88             :  * \p Faces of type Quad4.
      89             :  *
      90             :  * \author Benjamin S. Kirk
      91             :  * \date 2002-2007
      92             :  * \brief The base class for all geometric element types.
      93             :  */
      94             : class Elem : public ReferenceCountedObject<Elem>,
      95             :              public DofObject
      96             : {
      97             : protected:
      98             : 
      99             :   /**
     100             :    * Constructor.  Creates an element with \p n_nodes nodes,
     101             :    * \p n_sides sides, \p n_children possible children, and
     102             :    * parent \p p.  The constructor allocates the memory necessary
     103             :    * to support this data.
     104             :    */
     105             :   Elem (const unsigned int n_nodes,
     106             :         const unsigned int n_sides,
     107             :         Elem * parent,
     108             :         Elem ** elemlinkdata,
     109             :         Node ** nodelinkdata);
     110             : 
     111             : public:
     112             : 
     113             :   /**
     114             :    * Elems are responsible for allocating and deleting space for
     115             :    * storing pointers to their children during refinement, so they
     116             :    * cannot currently be (default) copy-constructed or copy-
     117             :    * assigned. We therefore explicitly delete these operations. In
     118             :    * addition, the DofObject base class currently has private copy
     119             :    * construction and assignment operators, so that prevents us from
     120             :    * copying Elems as well.
     121             :    */
     122             :   Elem (Elem &&) = delete;
     123             :   Elem (const Elem &) = delete;
     124             :   Elem & operator= (const Elem &) = delete;
     125             :   Elem & operator= (Elem &&) = delete;
     126             : 
     127             :   /**
     128             :    * Destructor.
     129             :    */
     130  3451074184 :   virtual ~Elem() = default;
     131             : 
     132             :   /**
     133             :    * \returns The \p Point associated with local \p Node \p i.
     134             :    */
     135             :   const Point & point (const unsigned int i) const;
     136             : 
     137             :   /**
     138             :    * \returns The \p Point associated with local \p Node \p i
     139             :    * as a writable reference.
     140             :    */
     141             :   Point & point (const unsigned int i);
     142             : 
     143             :   /**
     144             :    * \returns The \p Point associated with local \p Node \p i,
     145             :    * in master element rather than physical coordinates.
     146             :    */
     147             :   virtual Point master_point (const unsigned int i) const = 0;
     148             : 
     149             :   /**
     150             :    * \returns The global id number of local \p Node \p i.
     151             :    */
     152             :   dof_id_type node_id (const unsigned int i) const;
     153             : 
     154             :   /**
     155             :    * \returns The local id number of global \p Node id \p i,
     156             :    * or \p invalid_uint if Node id \p i is not local.
     157             :    */
     158             :   unsigned int local_node (const dof_id_type i) const;
     159             : 
     160             :   /**
     161             :    * \returns The local index for the \p Node pointer \p node_ptr,
     162             :    * or \p invalid_uint if \p node_ptr is not a local node.
     163             :    */
     164             :   unsigned int get_node_index (const Node * node_ptr) const;
     165             : 
     166             :   /**
     167             :    * \returns A pointer to an array of local node pointers.
     168             :    */
     169             :   const Node * const * get_nodes () const;
     170             : 
     171             :   /**
     172             :    * \returns A const pointer to local \p Node \p i.
     173             :    */
     174             :   const Node * node_ptr (const unsigned int i) const;
     175             : 
     176             :   /**
     177             :    * \returns A non-const pointer to local \p Node \p i.
     178             :    */
     179             :   Node * node_ptr (const unsigned int i);
     180             : 
     181             :   /**
     182             :    * \returns A const reference to local \p Node \p i.
     183             :    */
     184             :   const Node & node_ref (const unsigned int i) const;
     185             : 
     186             :   /**
     187             :    * \returns A writable reference to local \p Node \p i.
     188             :    */
     189             :   Node & node_ref (const unsigned int i);
     190             : 
     191             : #ifdef LIBMESH_ENABLE_DEPRECATED
     192             :   /**
     193             :    * \returns The pointer to the \p Node with local number \p i as a
     194             :    * writable reference.
     195             :    *
     196             :    * \deprecated This setter cannot update the multiple node pointers
     197             :    * used in a general polyhedron; use the \p set_node overload that
     198             :    * takes an argument.
     199             :    */
     200             :   virtual Node * & set_node (const unsigned int i);
     201             : #endif // LIBMESH_ENABLE_DEPRECATED
     202             : 
     203             :   /**
     204             :    * Sets local \p Node \p i to refer to \p node.
     205             :    */
     206             :   virtual void set_node (const unsigned int i,
     207             :                          Node * node);
     208             : 
     209             :   /**
     210             :    * Nested classes for use iterating over all nodes of an element.
     211             :    */
     212             :   class NodeRefIter;
     213             :   class ConstNodeRefIter;
     214             : 
     215             :   /**
     216             :    * Returns a range with all nodes of an element, usable in
     217             :    * range-based for loops.  The exact type of the return value here
     218             :    * may be subject to change in future libMesh releases, but the
     219             :    * iterators will always dereference to produce a reference to a
     220             :    * Node.
     221             :    */
     222             :   SimpleRange<NodeRefIter> node_ref_range();
     223             : 
     224             :   SimpleRange<ConstNodeRefIter> node_ref_range() const;
     225             : 
     226             :   /**
     227             :    * \returns The subdomain that this element belongs to.
     228             :    */
     229             :   subdomain_id_type subdomain_id () const;
     230             : 
     231             :   /**
     232             :    * \returns The subdomain that this element belongs to as a
     233             :    * writable reference.
     234             :    */
     235             :   subdomain_id_type & subdomain_id ();
     236             : 
     237             :   /**
     238             :    * A static integral constant representing an invalid subdomain id.
     239             :    * See also DofObject::{invalid_id, invalid_unique_id, invalid_processor_id}.
     240             :    *
     241             :    * \note We don't use the static_cast(-1) trick here since
     242             :    * \p subdomain_id_type is sometimes a *signed* integer for
     243             :    * compatibility reasons (see libmesh/id_types.h).
     244             :    */
     245             :   static constexpr subdomain_id_type invalid_subdomain_id
     246             :     = std::numeric_limits<subdomain_id_type>::max();
     247             : 
     248             :   /**
     249             :    * \returns true iff this element type can vary in topology (e.g.
     250             :    * have different numbers of sides and/or nodes) at runtime.  For
     251             :    * such general polygons or polyhedra, APIs which assume a fixed
     252             :    * topology are not safe to use.
     253             :    */
     254   291556478 :   virtual bool runtime_topology() const { return false; }
     255             : 
     256             :   /**
     257             :    * \returns A pointer to the "reference element" associated
     258             :    * with this element.  The reference element is the image of this
     259             :    * element in reference parametric space. Importantly, it is *not*
     260             :    * an actual element in the mesh, but rather a Singleton-type
     261             :    * object, so for example all \p Quad4 elements share the same
     262             :    * \p reference_elem().
     263             :    *
     264             :    * If the element is of a type that can admit multiple topologies,
     265             :    * such as a Polygon subtype, then there is no reference element;
     266             :    * for such types this method should not be used.
     267             :    */
     268             :   const Elem * reference_elem () const;
     269             : 
     270             :   /**
     271             :    * \returns An id associated with the \p s side of this element.
     272             :    * The id is not necessarily unique, but should be close.
     273             :    */
     274             :   virtual dof_id_type key (const unsigned int s) const = 0;
     275             : 
     276             :   /**
     277             :    * \returns An id associated with the \p s side of this element, as
     278             :    * defined solely by element vertices.  The id is not necessarily
     279             :    * unique, but should be close.  This is particularly useful in the
     280             :    * \p MeshBase::find_neighbors() routine.
     281             :    */
     282             :   virtual dof_id_type low_order_key (const unsigned int s) const = 0;
     283             : 
     284             :   /**
     285             :    * \returns An id associated with the global node ids of this
     286             :    * element.  The id is not necessarily unique, but should be
     287             :    * close. Uses the same hash as the key(s) function, so for example
     288             :    * if "tri3" is side 0 of "tet4", then tri3->key()==tet4->key(0).
     289             :    */
     290             :   virtual dof_id_type key () const;
     291             : 
     292             :   /**
     293             :    * \returns \p true if two elements are equivalent, \p false
     294             :    * otherwise.  This is true if the elements are connected to
     295             :    * identical global nodes, regardless of how those nodes might be
     296             :    * numbered local to the elements.
     297             :    */
     298             :   bool operator == (const Elem & rhs) const;
     299             : 
     300             :   /**
     301             :    * \returns \p false if two elements are equivalent, \p true
     302             :    * otherwise.
     303             :    */
     304             :   bool operator != (const Elem & rhs) const;
     305             : 
     306             :   /**
     307             :    * \returns \p true if two elements have equal topologies, false
     308             :    * otherwise.
     309             :    * This is true if the elements connect to nodes of the same id in
     310             :    * the same order, and neighbors of the same id on each side, the
     311             :    * same id on any parent and/or interior_parent link, etc.
     312             :    */
     313             :   bool topologically_equal (const Elem & rhs) const;
     314             : 
     315             :   /**
     316             :    * \returns A const pointer to the \f$ i^{th} \f$ neighbor of this
     317             :    * element, or \p nullptr if \p MeshBase::find_neighbors() has not been
     318             :    * called.
     319             :    *
     320             :    * \note If \p MeshBase::find_neighbors() has been called and this
     321             :    * function still returns \p nullptr, then the side is on a boundary of
     322             :    * the domain.
     323             :    */
     324             :   const Elem * neighbor_ptr (unsigned int i) const;
     325             : 
     326             :   /**
     327             :    * \returns A non-const pointer to the \f$ i^{th} \f$ neighbor of this element.
     328             :    */
     329             :   Elem * neighbor_ptr (unsigned int i);
     330             : 
     331             :   /**
     332             :    * Nested "classes" for use iterating over all neighbors of an element.
     333             :    */
     334             :   typedef Elem * const *       NeighborPtrIter;
     335             :   typedef const Elem * const * ConstNeighborPtrIter;
     336             : 
     337             :   /**
     338             :    * Returns a range with all neighbors of an element, usable in
     339             :    * range-based for loops.  The exact type of the return value here
     340             :    * may be subject to change in future libMesh releases, but the
     341             :    * iterators will always dereference to produce a pointer to a
     342             :    * neighbor element (or a null pointer, for sides which have no
     343             :    * neighbors).
     344             :    */
     345             :   SimpleRange<NeighborPtrIter> neighbor_ptr_range();
     346             : 
     347             :   SimpleRange<ConstNeighborPtrIter> neighbor_ptr_range() const;
     348             : 
     349             : #ifdef LIBMESH_ENABLE_PERIODIC
     350             :   /**
     351             :    * \returns A pointer to the \f$ i^{th} \f$ neighbor of this element
     352             :    * for interior elements.  If an element is on a periodic
     353             :    * boundary, it will return a corresponding element on the opposite
     354             :    * side.
     355             :    */
     356             :   const Elem * topological_neighbor (const unsigned int i,
     357             :                                      const MeshBase & mesh,
     358             :                                      const PointLocatorBase & point_locator,
     359             :                                      const PeriodicBoundaries * pb) const;
     360             : 
     361             :   /**
     362             :    * \returns A writable pointer to the \f$ i^{th} \f$ neighbor of
     363             :    * this element for interior elements.  If an element is on a
     364             :    * periodic boundary, it will return a corresponding element on the
     365             :    * opposite side.
     366             :    */
     367             :   Elem * topological_neighbor (const unsigned int i,
     368             :                                MeshBase & mesh,
     369             :                                const PointLocatorBase & point_locator,
     370             :                                const PeriodicBoundaries * pb);
     371             : 
     372             :   /**
     373             :    * \returns \p true if the element \p elem in question is a neighbor or
     374             :    * topological neighbor of this element, \p false otherwise.
     375             :    */
     376             :   bool has_topological_neighbor (const Elem * elem,
     377             :                                  const MeshBase & mesh,
     378             :                                  const PointLocatorBase & point_locator,
     379             :                                  const PeriodicBoundaries * pb) const;
     380             : #endif
     381             : 
     382             :   /**
     383             :    * Assigns \p n as the \f$ i^{th} \f$ neighbor.
     384             :    */
     385             :   void set_neighbor (const unsigned int i, Elem * n);
     386             : 
     387             :   /**
     388             :    * \returns \p true if the element \p elem in question is a neighbor
     389             :    * of this element, \p false otherwise.
     390             :    */
     391             :   bool has_neighbor (const Elem * elem) const;
     392             : 
     393             :   /**
     394             :    * \returns If \p elem is a neighbor of a child of this element, a
     395             :    * pointer to that child, otherwise \p nullptr.
     396             :    */
     397             :   Elem * child_neighbor (Elem * elem);
     398             : 
     399             :   /**
     400             :    * \returns If \p elem is a neighbor of a child of this element, a
     401             :    * pointer to that child, otherwise \p nullptr.
     402             :    */
     403             :   const Elem * child_neighbor (const Elem * elem) const;
     404             : 
     405             :   /**
     406             :    * \returns \p true if this element has a side coincident
     407             :    * with a boundary (indicated by a \p nullptr neighbor), \p false
     408             :    * otherwise.
     409             :    */
     410             :   bool on_boundary () const;
     411             : 
     412             :   /**
     413             :    * \returns \p true if this element is "semilocal" to the calling
     414             :    * processor, which must specify its rank.
     415             :    *
     416             :    * This method is discouraged, as it uses the *old* definition of
     417             :    * semilocal (elements which are not local but which are point
     418             :    * neighbors of something local) rather than any of the new
     419             :    * definitions discussed in ghosting_functor.h
     420             :    */
     421             :   bool is_semilocal (const processor_id_type my_pid) const;
     422             : 
     423             :   /**
     424             :    * This function tells you which neighbor \p e is.
     425             :    * I.e. if s = a->which_neighbor_am_i(e); then
     426             :    * a->neighbor(s) will be an ancestor of e.
     427             :    */
     428             :   unsigned int which_neighbor_am_i(const Elem * e) const;
     429             : 
     430             :   /**
     431             :    * This function tells you which side the boundary element \p e is.
     432             :    * I.e. if e = a->build_side_ptr(s) or e = a->side_ptr(s); then
     433             :    * a->which_side_am_i(e) will be s.
     434             :    *
     435             :    * \note An \e exact floating point comparison of the nodal
     436             :    * positions of \p e is made with the nodal positions of \p this in
     437             :    * order to perform this test. The idea is that the test will return
     438             :    * a valid side id if \p e either directly shares Node pointers with
     439             :    * \p this, or was created by exactly copying some of the nodes of
     440             :    * \p this (e.g. through BoundaryMesh::sync()). In these
     441             :    * circumstances, non-fuzzy floating point equality is expected.
     442             :    *
     443             :    * \returns The side of \p this the element which \p e is, otherwise
     444             :    * \p invalid_uint.
     445             :    */
     446             :   unsigned int which_side_am_i(const Elem * e) const;
     447             : 
     448             :   /**
     449             :    * \returns The local node id for node \p side_node on side \p side of
     450             :    * this Elem. Simply relies on the \p side_nodes_map for each of the
     451             :    * derived types. For example,
     452             :    * Tri3::local_side_node(0, 0) -> 0
     453             :    * Tri3::local_side_node(0, 1) -> 1
     454             :    * Tri3::local_side_node(1, 0) -> 1
     455             :    * Tri3::local_side_node(1, 1) -> 2
     456             :    * etc...
     457             :    */
     458             :   virtual unsigned int local_side_node(unsigned int side,
     459             :                                        unsigned int side_node) const = 0;
     460             : 
     461             :   /**
     462             :    * Similar to Elem::local_side_node(), but instead of a side id, takes
     463             :    * an edge id and a node id on that edge and returns a local node number
     464             :    * for the Elem. The implementation relies on the "edge_nodes_map" tables
     465             :    * for 3D elements. For 2D elements, calls local_side_node(). Throws an
     466             :    * error if called on 1D elements.
     467             :    */
     468             :   virtual unsigned int local_edge_node(unsigned int edge,
     469             :                                        unsigned int edge_node) const = 0;
     470             : 
     471             :   /**
     472             :    * \returns \p true if a vertex of \p e is contained
     473             :    * in this element.  If \p mesh_connection is true, looks
     474             :    * specifically for containment possibilities of an element \p e
     475             :    * that is connected to \p this via membership in the same manifold
     476             :    * of the same mesh.
     477             :    */
     478             :   bool contains_vertex_of(const Elem * e, bool mesh_connection=false) const;
     479             : 
     480             :   /**
     481             :    * \returns \p true if an edge of \p e is contained in
     482             :    * this element.  (Internally, this is done by checking whether at
     483             :    * least two vertices of \p e are contained in this element).
     484             :    */
     485             :   bool contains_edge_of(const Elem * e) const;
     486             : 
     487             :   /**
     488             :    * This function finds all active elements (including this one)
     489             :    * which are in the same manifold as this element and which touch
     490             :    * the current active element at the specified point, which should
     491             :    * be a point in the current element.
     492             :    *
     493             :    * Elements which are not "in the same manifold" (e.g. the
     494             :    * interior_parent of a boundary element) will not be found with
     495             :    * this method.
     496             :    *
     497             :    * Elements which overlap the specified point but which are only
     498             :    * connected to the current element via elements which do not
     499             :    * overlap that point (e.g. in a folded or tangled mesh) are not
     500             :    * considered to "touch" the current element and will not be found
     501             :    * with this method.
     502             :    */
     503             :   void find_point_neighbors(const Point & p,
     504             :                             std::set<const Elem *> & neighbor_set) const;
     505             : 
     506             :   /**
     507             :    * This function finds all active elements (including this one) in
     508             :    * the same manifold as this element which touch this active element
     509             :    * at any point.
     510             :    */
     511             :   void find_point_neighbors(std::set<const Elem *> & neighbor_set) const;
     512             : 
     513             :   /**
     514             :    * This function finds all active elements (including this one) in
     515             :    * the same manifold as start_elem (which must be active and must
     516             :    * touch this element) which touch this element at any point.
     517             :    */
     518             :   void find_point_neighbors(std::set<const Elem *> & neighbor_set,
     519             :                             const Elem * start_elem) const;
     520             : 
     521             :   /**
     522             :    * Non-const version of function above. Fills a set of non-const Elem pointers.
     523             :    */
     524             :   void find_point_neighbors(std::set<Elem *> & neighbor_set,
     525             :                             Elem * start_elem);
     526             : 
     527             :   /**
     528             :    * This function finds all active elements in the same manifold as
     529             :    * this element which touch the current active element along the
     530             :    * whole edge defined by the two points \p p1 and \p p2.
     531             :    */
     532             :   void find_edge_neighbors(const Point & p1,
     533             :                            const Point & p2,
     534             :                            std::set<const Elem *> & neighbor_set) const;
     535             : 
     536             :   /**
     537             :    * This function finds all active elements in the same manifold as
     538             :    * this element which touch the current active element along any
     539             :    * edge (more precisely, at at least two points).
     540             :    *
     541             :    * In this case, elements are included even if they do not touch a
     542             :    * *whole* edge of this element.
     543             :    */
     544             :   void find_edge_neighbors(std::set<const Elem *> & neighbor_set) const;
     545             : 
     546             :   /**
     547             :    * This function finds all active elements (*not* including this
     548             :    * one) in the parent manifold of this element whose intersection
     549             :    * with this element has non-zero measure.
     550             :    */
     551             :   void find_interior_neighbors(std::set<const Elem *> & neighbor_set) const;
     552             : 
     553             :   /**
     554             :    * Non-const version of function above that fills up a vector of
     555             :    * non-const Elem pointers instead.
     556             :    */
     557             :   void find_interior_neighbors(std::set<Elem *> & neighbor_set);
     558             : 
     559             :   /**
     560             :    * Resets this element's neighbors' appropriate neighbor pointers
     561             :    * and its parent's and children's appropriate pointers
     562             :    * to point to null instead of to this.
     563             :    *
     564             :    * To be used before an element is deleted from a mesh.
     565             :    */
     566             :   void remove_links_to_me ();
     567             : 
     568             :   /**
     569             :    * Resets this element's neighbors' appropriate neighbor pointers
     570             :    * and its parent's and children's appropriate pointers
     571             :    * to point to the global remote_elem instead of this.
     572             :    * Used by the library before an element becomes remote on the
     573             :    * local processor.
     574             :    */
     575             :   void make_links_to_me_remote ();
     576             : 
     577             :   /**
     578             :    * Resets the \p neighbor_side pointers of our nth neighbor (and
     579             :    * its descendants, if appropriate) to point to this Elem instead of
     580             :    * to the global remote_elem.  Used by the library when a formerly
     581             :    * remote element is being added to the local processor.
     582             :    */
     583             :   void make_links_to_me_local (unsigned int n, unsigned int neighbor_side);
     584             : 
     585             :   /**
     586             :    * \returns \p true if this element is remote, false otherwise.
     587             :    *
     588             :    * A remote element (see \p RemoteElem) is a syntactic convenience --
     589             :    * it is a placeholder for an element which exists on some other
     590             :    * processor.  Local elements are required to have valid neighbors,
     591             :    * and these ghost elements may have remote neighbors for data
     592             :    * structure consistency.  The use of remote elements helps ensure
     593             :    * that any element we may access has a \p nullptr neighbor only if it
     594             :    * lies on the physical boundary of the domain.
     595             :    */
     596     7280126 :   virtual bool is_remote () const
     597     7280126 :   { return false; }
     598             : 
     599             :   /**
     600             :    * \returns The connectivity for this element in a specific
     601             :    * format, which is specified by the IOPackage tag.
     602             :    */
     603             :   virtual void connectivity(const unsigned int sc,
     604             :                             const IOPackage iop,
     605             :                             std::vector<dof_id_type> & conn) const = 0;
     606             : 
     607             :   /**
     608             :    * Writes the element connectivity for various IO packages
     609             :    * to the passed ostream "out".  Not virtual, since it is
     610             :    * implemented in the base class.
     611             :    */
     612             :   void write_connectivity (std::ostream & out,
     613             :                            const IOPackage iop) const;
     614             : 
     615             :   /**
     616             :    * \returns The type of element that has been derived from this
     617             :    * base class.
     618             :    */
     619             :   virtual ElemType type () const = 0;
     620             : 
     621             :   /**
     622             :    * This array maps the integer representation of the \p ElemType enum
     623             :    * to the geometric dimension of the element.
     624             :    *
     625             :    * This is currently usable even for complicated subclasses with
     626             :    * runtime-varying topology.
     627             :    */
     628             :   static const unsigned int type_to_dim_map[INVALID_ELEM];
     629             : 
     630             :   /**
     631             :    * \returns The dimensionality of the object.
     632             :    */
     633             :   virtual unsigned short dim () const = 0;
     634             : 
     635             :   /**
     636             :    * This array maps the integer representation of the \p ElemType enum
     637             :    * to the number of nodes in the element.
     638             :    *
     639             :    * This is only usable for simple types for which the node number
     640             :    * is fixed; for more general types like Polygon subclasses an actual
     641             :    * instantiated Elem must be queried.
     642             :    */
     643             :   static const unsigned int type_to_n_nodes_map[INVALID_ELEM];
     644             : 
     645             :   /**
     646             :    * \returns The number of nodes this element contains.
     647             :    */
     648             :   virtual unsigned int n_nodes () const = 0;
     649             : 
     650             :   /**
     651             :    * The maximum number of nodes *any* element can contain.
     652             :    * This is useful for replacing heap vectors with stack arrays.
     653             :    */
     654             :   static const unsigned int max_n_nodes = 27;
     655             : 
     656             :   /**
     657             :    * \returns An integer range from 0 up to (but not including)
     658             :    * the number of nodes this element contains.
     659             :    */
     660             :   IntRange<unsigned short> node_index_range () const;
     661             : 
     662             :   /**
     663             :    * \returns The number of nodes the given child of this element
     664             :    * contains.  Except in odd cases like pyramid refinement this will
     665             :    * be the same as the number of nodes in the parent element.
     666             :    */
     667    21575143 :   virtual unsigned int n_nodes_in_child (unsigned int /*c*/) const
     668    21575143 :   { return this->n_nodes(); }
     669             : 
     670             :   /**
     671             :    * This array maps the integer representation of the \p ElemType enum
     672             :    * to the number of sides on the element.
     673             :    *
     674             :    * This is only usable for simple types for which the node number
     675             :    * is fixed; for more general types like Polygon subclasses an actual
     676             :    * instantiated Elem must be queried.
     677             :    */
     678             :   static const unsigned int type_to_n_sides_map[INVALID_ELEM];
     679             : 
     680             :   /**
     681             :    * \returns The number of sides the element that has been derived
     682             :    * from this class has. In 2D the number of sides is the number
     683             :    * of edges, in 3D the number of sides is the number of faces.
     684             :    */
     685             :   virtual unsigned int n_sides () const = 0;
     686             : 
     687             :   /**
     688             :    * \returns The type of element for side \p s.
     689             :    */
     690             :   virtual ElemType side_type (const unsigned int s) const = 0;
     691             : 
     692             :   /**
     693             :    * \returns the normal (outwards-facing) of the side of the element at the vertex-average of the side
     694             :    * @param s the side of interest
     695             :    */
     696             :   virtual Point side_vertex_average_normal(const unsigned int s) const;
     697             : 
     698             :   /**
     699             :    * \returns An integer range from 0 up to (but not including)
     700             :    * the number of sides this element has.
     701             :    */
     702             :   IntRange<unsigned short> side_index_range () const;
     703             : 
     704             :   /**
     705             :    * \returns The number of neighbors the element that has been derived
     706             :    * from this class has.
     707             :    *
     708             :    * Only face (or edge in 2D) neighbors are stored, so this method
     709             :    * returns n_sides().  At one point we intended to allow derived
     710             :    * classes to override this, but too much current libMesh code
     711             :    * assumes n_neighbors==n_sides.
     712             :    */
     713    84694137 :   unsigned int n_neighbors () const
     714   161549458 :   { return this->n_sides(); }
     715             : 
     716             :   /**
     717             :    * \returns The number of vertices the element that has been derived
     718             :    * from this class has.
     719             :    */
     720             :   virtual unsigned int n_vertices () const = 0;
     721             : 
     722             :   /**
     723             :    * \returns The number of edges the element that has been derived
     724             :    * from this class has.
     725             :    */
     726             :   virtual unsigned int n_edges () const = 0;
     727             : 
     728             :   /**
     729             :    * \returns An integer range from 0 up to (but not including)
     730             :    * the number of edges this element has.
     731             :    */
     732             :   IntRange<unsigned short> edge_index_range () const;
     733             : 
     734             :   /**
     735             :    * This array maps the integer representation of the \p ElemType enum
     736             :    * to the number of edges on the element.
     737             :    *
     738             :    * This is only usable for simple types for which the node number
     739             :    * is fixed; for more general types like Polygon subclasses an actual
     740             :    * instantiated Elem must be queried.
     741             :    */
     742             :   static const unsigned int type_to_n_edges_map[INVALID_ELEM];
     743             : 
     744             :   /**
     745             :    * \returns The number of faces the element that has been derived
     746             :    * from this class has.
     747             :    */
     748             :   virtual unsigned int n_faces () const = 0;
     749             : 
     750             :   /**
     751             :    * \returns An integer range from 0 up to (but not including)
     752             :    * the number of faces this element has.
     753             :    */
     754             :   IntRange<unsigned short> face_index_range () const;
     755             : 
     756             :   /**
     757             :    * \returns The number of children the element that has been derived
     758             :    * from this class may have.
     759             :    */
     760             :   virtual unsigned int n_children () const = 0;
     761             : 
     762             :   /**
     763             :    * \returns \p true if the specified (local) node number is a vertex node.
     764             :    */
     765             :   virtual bool is_vertex(const unsigned int i) const = 0;
     766             : 
     767             :   /**
     768             :    * \returns \p true if the specified child has a vertex at the
     769             :    * specified (child-local) node number.
     770             :    * Except in odd cases like pyramid refinement the child will have
     771             :    * the same local structure as the parent element.
     772             :    */
     773      177423 :   virtual bool is_vertex_on_child (unsigned int /*c*/,
     774             :                                    unsigned int n) const
     775      177423 :   { return this->is_vertex(n); }
     776             : 
     777             :   /**
     778             :    * \returns \p true if this element has a vertex at the specified
     779             :    * (child-local) node number \p n of the specified child \p c.
     780             :    */
     781             :   virtual bool is_vertex_on_parent(unsigned int c,
     782             :                                    unsigned int n) const;
     783             : 
     784             :   /**
     785             :    * \returns \p true if the specified (local) node number is an edge node.
     786             :    * For 1D elements, is_edge() is equivalent to is_internal().
     787             :    */
     788             :   virtual bool is_edge(const unsigned int i) const = 0;
     789             : 
     790             :   /**
     791             :    * \returns \p true if the specified (local) node number is a face node.
     792             :    * For 2D elements, is_face() is equivalent to is_internal().
     793             :    * For 1D elements, is_face() == false.
     794             :    */
     795             :   virtual bool is_face(const unsigned int i) const = 0;
     796             : 
     797             :   /**
     798             :    * \returns \p true if the specified (local) node number is an internal node.
     799             :    */
     800             :   bool is_internal(const unsigned int i) const;
     801             : 
     802             :   /**
     803             :    * \returns \p true if the specified (local) node number is on the
     804             :    * specified side.
     805             :    */
     806             :   virtual bool is_node_on_side(const unsigned int n,
     807             :                                const unsigned int s) const = 0;
     808             : 
     809             :   /**
     810             :    * \returns the (local) node numbers on the specified side
     811             :    */
     812             :   virtual std::vector<unsigned int> nodes_on_side(const unsigned int /*s*/) const = 0;
     813             : 
     814             :   /**
     815             :    * \returns the (local) node numbers on the specified edge
     816             :    */
     817             :   virtual std::vector<unsigned int> nodes_on_edge(const unsigned int /*e*/) const = 0;
     818             : 
     819             :   /**
     820             :    * \returns the (local) side numbers that touch the specified edge
     821             :    */
     822             :   virtual std::vector<unsigned int> sides_on_edge(const unsigned int /*e*/) const = 0;
     823             : 
     824             :   /**
     825             :    * \returns the (local) edge numbers that touch the specified node
     826             :    */
     827             :   virtual std::vector<unsigned int> edges_adjacent_to_node(const unsigned int /*n*/) const = 0;
     828             : 
     829             :   /**
     830             :    * \returns \p true if the specified (local) node number is on the
     831             :    * specified edge.
     832             :    */
     833             :   virtual bool is_node_on_edge(const unsigned int n,
     834             :                                const unsigned int e) const = 0;
     835             : 
     836             :   /**
     837             :    * \returns \p true if the specified edge is on the specified side.
     838             :    */
     839             :   virtual bool is_edge_on_side(const unsigned int e,
     840             :                                const unsigned int s) const = 0;
     841             : 
     842             :   /**
     843             :    * \returns The side number opposite to \p s (for a tensor product
     844             :    * element), or throws an error otherwise.
     845             :    */
     846             :   virtual unsigned int opposite_side(const unsigned int s) const;
     847             : 
     848             :   /**
     849             :    * \returns The local node number for the node opposite to node n
     850             :    * on side \p opposite_side(s) (for a tensor product element), or
     851             :    * throws an error otherwise.
     852             :    */
     853             :   virtual unsigned int opposite_node(const unsigned int n,
     854             :                                      const unsigned int s) const;
     855             : 
     856             :   /**
     857             :    * \returns The number of sub-elements this element may be broken
     858             :    * down into for visualization purposes.  For example, 1 for a
     859             :    * linear triangle, 4 for a quadratic (6-noded) triangle, etc...
     860             :    */
     861             :   virtual unsigned int n_sub_elem () const = 0;
     862             : 
     863             :   /**
     864             :    * \returns A temporary element coincident with side \p i.
     865             :    *
     866             :    * This method returns the _minimum_ element necessary to uniquely
     867             :    * identify the side.  For example, the side of a hexahedron is
     868             :    * always returned as a 4-noded quadrilateral, regardless of what
     869             :    * type of hex you are dealing with.  Important data like subdomain
     870             :    * id, p level, or mapping type may be omitted from the temporary
     871             :    * element.  If you want a first-class full-ordered face (i.e. a
     872             :    * 9-noded quad face for a 27-noded hexahedron), use the
     873             :    * build_side_ptr method.
     874             :    *
     875             :    * \note The const version of this function is non-virtual; it
     876             :    * simply calls the virtual non-const version and const_casts the
     877             :    * return type.
     878             :    */
     879             :   virtual std::unique_ptr<Elem> side_ptr (unsigned int i) = 0;
     880             :   std::unique_ptr<const Elem> side_ptr (unsigned int i) const;
     881             : 
     882             :   /**
     883             :    * Resets the loose element \p side, which may currently point to a
     884             :    * different side than \p i or even a different element than \p
     885             :    * this, to point to side \p i on \p this.  If \p side is currently
     886             :    * an element of the wrong type, it will be freed and a new element
     887             :    * allocated; otherwise no memory allocation will occur.
     888             :    *
     889             :    * This will cause \p side to be a minimum-ordered element, even if
     890             :    * it is handed a higher-ordered element that must be replaced.
     891             :    *
     892             :    * The const version of this function is non-virtual; it simply
     893             :    * calls the virtual non-const version and const_casts the return
     894             :    * type.
     895             :    */
     896             :   virtual void side_ptr (std::unique_ptr<Elem> & side, const unsigned int i) = 0;
     897             :   void side_ptr (std::unique_ptr<const Elem> & side, const unsigned int i) const;
     898             : 
     899             :   /**
     900             :    * \returns An temporary element coincident with side \p i wrapped
     901             :    * in a smart pointer.
     902             :    *
     903             :    * The element returned is full-ordered and full-featured, in
     904             :    * contrast to the side method.  For example, calling
     905             :    * build_side_ptr(0) on a 20-noded hex in subdomain 5 will build a
     906             :    * 8-noded quadrilateral coincident with face 0, assign it subdomain
     907             :    * id 5, and pass back the pointer.
     908             :    *
     909             :    * The side element's id() is undefined; it is a temporary element
     910             :    * not added to any mesh.
     911             :    *
     912             :    * A \p std::unique_ptr<Elem> is returned to prevent a memory leak.
     913             :    * This way the user need not remember to delete the object.
     914             :    *
     915             :    * The const version of this function is non-virtual; it simply
     916             :    * calls the virtual non-const version and const_casts the return
     917             :    * type.
     918             :    */
     919             :   virtual std::unique_ptr<Elem> build_side_ptr (const unsigned int i) = 0;
     920             :   std::unique_ptr<const Elem> build_side_ptr (const unsigned int i) const;
     921             : 
     922             : #ifdef LIBMESH_ENABLE_DEPRECATED
     923             :   /*
     924             :    * Older versions of libMesh supported a "proxy" option here.
     925             :    */
     926           0 :   virtual std::unique_ptr<Elem> build_side_ptr (const unsigned int i, bool proxy)
     927           0 :   { if (proxy) libmesh_error(); libmesh_deprecated(); return this->build_side_ptr(i); }
     928             : 
     929             :   std::unique_ptr<const Elem> build_side_ptr (const unsigned int i, bool proxy) const
     930             :   { if (proxy) libmesh_error(); libmesh_deprecated(); return this->build_side_ptr(i); }
     931             : #endif
     932             : 
     933             :   /**
     934             :    * Resets the loose element \p side, which may currently point to a
     935             :    * different side than \p i or even a different element than \p
     936             :    * this, to point to side \p i on \p this.  If \p side is currently
     937             :    * an element of the wrong type, it will be freed and a new element
     938             :    * allocated; otherwise no memory allocation will occur.
     939             :    *
     940             :    * This will cause \p side to be a full-ordered element, even if it
     941             :    * is handed a lower-ordered element that must be replaced.
     942             :    *
     943             :    * The const version of this function is non-virtual; it simply
     944             :    * calls the virtual non-const version and const_casts the return
     945             :    * type.
     946             :    */
     947             :   virtual void build_side_ptr (std::unique_ptr<Elem> & side, const unsigned int i) = 0;
     948             :   void build_side_ptr (std::unique_ptr<const Elem> & side, const unsigned int i) const;
     949             : 
     950             :   /**
     951             :    * \returns An element coincident with edge \p i wrapped in a smart pointer.
     952             :    *
     953             :    * The element returned is full-ordered.  For example, calling
     954             :    * build_edge_ptr(0) on a 20-noded hex will build a 3-noded edge
     955             :    * coincident with edge 0 and pass back the pointer.  A \p
     956             :    * std::unique_ptr<Elem> is returned to prevent a memory leak.  This way
     957             :    * the user need not remember to delete the object.
     958             :    *
     959             :    * The const version of this function is non-virtual; it simply
     960             :    * calls the virtual non-const version and const_casts the return
     961             :    * type.
     962             :    */
     963             :   virtual std::unique_ptr<Elem> build_edge_ptr (const unsigned int i) = 0;
     964             :   std::unique_ptr<const Elem> build_edge_ptr (const unsigned int i) const;
     965             : 
     966             :   /**
     967             :    * Resets the loose element \p edge, which may currently point to a
     968             :    * different edge than \p i or even a different element than \p
     969             :    * this, to point to edge \p i on \p this.  If \p edge is currently
     970             :    * an element of the wrong type, it will be freed and a new element
     971             :    * allocated; otherwise no memory allocation will occur.
     972             :    *
     973             :    * This will cause \p edge to be a full-ordered element, even if it
     974             :    * is handed a lower-ordered element that must be replaced.
     975             :    *
     976             :    * The const version of this function is non-virtual; it simply
     977             :    * calls the virtual non-const version and const_casts the return
     978             :    * type.
     979             :    */
     980             :   virtual void build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i) = 0;
     981             :   void build_edge_ptr (std::unique_ptr<const Elem> & edge, const unsigned int i) const;
     982             : 
     983             :   /**
     984             :    * This array maps the integer representation of the \p ElemType enum
     985             :    * to the default approximation order of elements of that type.
     986             :    *
     987             :    * This is currently usable even for complicated subclasses with
     988             :    * runtime-varying topology.
     989             :    */
     990             :   static const Order type_to_default_order_map[INVALID_ELEM];
     991             : 
     992             :   /**
     993             :    * \returns The default approximation order for this element.  This
     994             :    * is the order that will be used to compute the map to the
     995             :    * reference element.
     996             :    */
     997             :   virtual Order default_order () const = 0;
     998             : 
     999             :   /**
    1000             :    * \returns The maximum supported approximation order for nodal
    1001             :    * (Lagrange or Rational Bezier-Bernstein) variables on this element
    1002             :    * type.  This is usually the same as the default order.
    1003             :    */
    1004    31856884 :   virtual Order supported_nodal_order() const { return default_order(); }
    1005             : 
    1006             :   /**
    1007             :    * \returns The default approximation order for side elements of
    1008             :    * this element type.  This may be lower for elements with 'bubble
    1009             :    * functions' in the Lagrange basis.
    1010             :    */
    1011      513235 :   virtual Order default_side_order () const { return default_order(); }
    1012             : 
    1013             :   /**
    1014             :    * \returns The "true" geometric centroid of the element, c=(cx, cy,
    1015             :    * cz), where:
    1016             :    *
    1017             :    * [cx]            [\int x dV]
    1018             :    * [cy] := (1/V) * [\int y dV]
    1019             :    * [cz]            [\int z dV]
    1020             :    *
    1021             :    * This method is virtual since some derived elements might want to
    1022             :    * use shortcuts to compute their centroid. For most element types,
    1023             :    * this method is more expensive than calling vertex_average(), so
    1024             :    * if you only need a point which is located "somewhere" in the
    1025             :    * interior of the element, consider calling vertex_average() instead.
    1026             :    */
    1027             :   virtual Point true_centroid () const;
    1028             : 
    1029             :   /**
    1030             :    * \returns A Point at the average of the elment's vertices.
    1031             :    *
    1032             :    * \note This used to be the base class centroid() implementation, but
    1033             :    * the centroid is only equal to the vertex average in some special cases.
    1034             :    * The centroid() implementation now returns the "true" centroid of the
    1035             :    * element (up to quadrature error).
    1036             :    */
    1037             :   Point vertex_average () const;
    1038             : 
    1039             :   /**
    1040             :    * \returns The "circumcenter of mass" (area-weighted average of
    1041             :    * triangulation circumcenters) of the element.
    1042             :    *
    1043             :    * Not implemented for infinite elements, not currently implemented
    1044             :    * for 3D elements, currently ignores curvature of element edges.
    1045             :    */
    1046           0 :   virtual Point quasicircumcenter () const
    1047           0 :   { libmesh_not_implemented(); }
    1048             : 
    1049             :   /**
    1050             :    * \returns The minimum vertex separation for the element.
    1051             :    */
    1052             :   virtual Real hmin () const;
    1053             : 
    1054             :   /**
    1055             :    * \returns The maximum vertex separation for the element.
    1056             :    */
    1057             :   virtual Real hmax () const;
    1058             : 
    1059             :   /**
    1060             :    * \returns The (length/area/volume) of the geometric element.
    1061             :    *
    1062             :    * If the element is twisted or inverted such that the mapping
    1063             :    * Jacobian is singular at any point, implementations of this method
    1064             :    * may return a "net" volume or may simply return NaN.
    1065             :    */
    1066             :   virtual Real volume () const;
    1067             : 
    1068             :   /**
    1069             :    * \returns A bounding box (not necessarily the minimal bounding box)
    1070             :    * containing the geometric element.
    1071             :    *
    1072             :    * The base class implementation determines a bounding box for the
    1073             :    * element *nodes*, which should be sufficient for first order
    1074             :    * finite elements.  Higher order geometric elements will need to
    1075             :    * override with an implementation which takes curved elements into
    1076             :    * account.
    1077             :    */
    1078             :   virtual BoundingBox loose_bounding_box () const;
    1079             : 
    1080             :   /**
    1081             :    * \returns A quantitative assessment of element quality based on
    1082             :    * the quality metric \p q specified by the user. Not all ElemQuality
    1083             :    * metrics are supported for all Elem types; consult the Elem::quality()
    1084             :    * overrides for specific Elem types to determine which quality metrics
    1085             :    * are supported. The ElemQuality metrics with generic support for all
    1086             :    * Elems with dimension > 1 are:
    1087             :    * .) EDGE_LENGTH_RATIO - ratio of maximum to minimum edge (in 2D,
    1088             :    *    side) length, where the min/max is taken over all Elem edges.
    1089             :    * .) MIN,MAX_ANGLE - The minimum (respectively maximum) angle
    1090             :    *    between all pairs of adjacent Elem edges, in degrees. Note
    1091             :    *    that, in 3D, these are *not* the dihedral angles (angle
    1092             :    *    between adjacent planar faces of the element), which we plan to
    1093             :    *    add support for in the future. In 2D, we compute the angle
    1094             :    *    between adjacent sides for this metric.
    1095             :    */
    1096             :   virtual Real quality (const ElemQuality q) const;
    1097             : 
    1098             :   /**
    1099             :    * \returns The suggested quality bounds for the Elem based on
    1100             :    * quality measure \p q.
    1101             :    *
    1102             :    * These are the values suggested by the CUBIT User's Manual.  Since
    1103             :    * this function can have no possible meaning for an abstract Elem,
    1104             :    * it is an error in the base class.
    1105             :    */
    1106           0 :   virtual std::pair<Real,Real> qual_bounds (const ElemQuality) const
    1107           0 :   { libmesh_not_implemented(); return std::make_pair(0.,0.); }
    1108             : 
    1109             :   /**
    1110             :    * \returns \p true if the physical point p is contained in this
    1111             :    * element, false otherwise.
    1112             :    *
    1113             :    * For linear elements, performs an initial tight bounding box check
    1114             :    * (as an optimization step) and (if that passes) then uses the
    1115             :    * user-defined tolerance "tol" in a call to inverse_map() to actually
    1116             :    * test if the point is in the element.  For quadratic elements, the
    1117             :    * bounding box optimization is skipped, and only the inverse_map()
    1118             :    * steps are performed.
    1119             :    *
    1120             :    * \note This routine should not be used to determine if a point
    1121             :    * is merely "nearby" an element to within some tolerance. For that,
    1122             :    * use Elem::close_to_point() instead.
    1123             :    */
    1124             :   virtual bool contains_point (const Point & p, Real tol=TOLERANCE) const;
    1125             : 
    1126             :   /**
    1127             :    * \returns \p true if the master-space point p is contained in the
    1128             :    * reference element corresponding to this element, false otherwise.
    1129             :    *
    1130             :    * Since we are doing floating point comparisons here the parameter
    1131             :    * \p eps can be specified to indicate a tolerance.  For example,
    1132             :    * \f$ x \le 1 \f$  becomes \f$ x \le 1 + \epsilon \f$.
    1133             :    */
    1134             :   virtual bool on_reference_element(const Point & p,
    1135             :                                     const Real eps = TOLERANCE) const = 0;
    1136             : 
    1137             :   /**
    1138             :    * \returns \p true if this element is "close" to the point p, where
    1139             :    * "close" is determined by the tolerance tol.
    1140             :    */
    1141             :   virtual bool close_to_point(const Point & p, Real tol) const;
    1142             : 
    1143             :   /**
    1144             :    * \returns \p true if edge \p i is positively oriented. An edge is
    1145             :    * positively oriented iff its first vertex (i.e. zeroth node) is
    1146             :    * lexicographically greater than its second vertex (i.e. first node).
    1147             :    */
    1148             :   bool positive_edge_orientation(const unsigned int i) const;
    1149             : 
    1150             :   /**
    1151             :    * \returns \p true if face \p i is positively oriented. A face is
    1152             :    * positively oriented iff the triangle defined by the lexicographically
    1153             :    * least vertex and its two adjacent vertices on the same face is
    1154             :    * positively oriented. Said triangle is positively oriented iff its
    1155             :    * vertices are an odd permutation of their lexicographic ordering.
    1156             :    */
    1157             :   bool positive_face_orientation(const unsigned int i) const;
    1158             : 
    1159             :   /**
    1160             :    * \returns \p true iff, for an edge \p e on side \p s, the node map for
    1161             :    * side \p s is such that the first vertex (i.e. zeroth node) of \p e is
    1162             :    * lower positioned than the second vertex (i.e. first node) of \p e.
    1163             :    */
    1164             :   bool relative_edge_face_order(const unsigned int e, const unsigned int s) const;
    1165             : 
    1166             :   /**
    1167             :    * A helper function for copying generic element data (mapping,
    1168             :    * subdomain, processor) from an element to a derived (child, side,
    1169             :    * edge) element.  Useful for forwards compatibility when new data
    1170             :    * is added.
    1171             :    */
    1172             :   void inherit_data_from(const Elem & src);
    1173             : 
    1174             : private:
    1175             :   /**
    1176             :    * Shared private implementation used by the contains_point()
    1177             :    * and close_to_point() routines.  The box_tol tolerance is
    1178             :    * used in the bounding box optimization, the map_tol tolerance is used
    1179             :    * in the calls to inverse_map() and on_reference_element().
    1180             :    */
    1181             :   bool point_test(const Point & p, Real box_tol, Real map_tol) const;
    1182             : 
    1183             : public:
    1184             :   /**
    1185             :    * \returns \p true if the element map is definitely affine (i.e. the same at
    1186             :    * every quadrature point) within numerical tolerances.
    1187             :    */
    1188           0 :   virtual bool has_affine_map () const { return false; }
    1189             : 
    1190             :   /**
    1191             :    * \returns \p true if the element map is invertible everywhere on
    1192             :    * the element, to within a user-specified tolerance. The tolerance
    1193             :    * is generally used in comparisons against zero, so it should be an
    1194             :    * absolute rather than a relative tolerance. Throws a
    1195             :    * libmesh_not_implemented() error unless specialized by derived
    1196             :    * classes.
    1197             :    */
    1198             :   virtual bool has_invertible_map(Real tol = TOLERANCE*TOLERANCE) const;
    1199             : 
    1200             :   /**
    1201             :    * \returns \p true if the Lagrange shape functions on this element
    1202             :    * are linear.
    1203             :    */
    1204           0 :   virtual bool is_linear () const { return false; }
    1205             : 
    1206             :   /**
    1207             :    * Prints relevant information about the element.
    1208             :    */
    1209             :   void print_info (std::ostream & os=libMesh::out) const;
    1210             : 
    1211             :   /**
    1212             :    * Prints relevant information about the element to a string.
    1213             :    */
    1214             :   std::string get_info () const;
    1215             : 
    1216             :   /**
    1217             :    * \returns \p true if the element is active (i.e. has no active
    1218             :    * descendants) or AMR is disabled, \p false otherwise.
    1219             :    *
    1220             :    * \note It suffices to check the first child only.
    1221             :    */
    1222             :   bool active () const;
    1223             : 
    1224             :   /**
    1225             :    * \returns \p true if the element is an ancestor (i.e. has an
    1226             :    * active child or ancestor child), \p false otherwise or when AMR
    1227             :    * is disabled.
    1228             :    */
    1229             :   bool ancestor () const;
    1230             : 
    1231             :   /**
    1232             :    * \returns \p true if the element is subactive (i.e. has no active
    1233             :    * descendants), \p false otherwise or if AMR is disabled.
    1234             :    */
    1235             :   bool subactive () const;
    1236             : 
    1237             :   /**
    1238             :    * \returns \p true if the element has any children (active or not),
    1239             :    * \p false otherwise, or if AMR is disabled.
    1240             :    */
    1241             :   bool has_children () const;
    1242             : 
    1243             :   /**
    1244             :    * \returns \p true if the element has any descendants other than
    1245             :    * its immediate children, \p false otherwise, or if AMR is disabled.
    1246             :    */
    1247             :   bool has_ancestor_children () const;
    1248             : 
    1249             :   /**
    1250             :    * \returns \p true if \p descendant is a child of \p this, or a
    1251             :    * child of a child of \p this, etc., \p false otherwise or if AMR
    1252             :    * is disabled.
    1253             :    */
    1254             :   bool is_ancestor_of(const Elem * descendant) const;
    1255             : 
    1256             :   /**
    1257             :    * \returns A const pointer to the element's parent, or \p nullptr if
    1258             :    * the element was not created via refinement.
    1259             :    */
    1260             :   const Elem * parent () const;
    1261             : 
    1262             :   /**
    1263             :    * \returns A pointer to the element's parent, or \p nullptr if
    1264             :    * the element was not created via refinement.
    1265             :    */
    1266             :   Elem * parent ();
    1267             : 
    1268             :   /**
    1269             :    * Sets the pointer to the element's parent.
    1270             :    * Dangerous! Only use this if you know what you are doing!
    1271             :    */
    1272             :   void set_parent (Elem * p);
    1273             : 
    1274             :   /**
    1275             :    * \returns A pointer to the element's top-most (i.e. level-0) parent.
    1276             :    *
    1277             :    * That is, \p this if this is a level-0 element, this element's parent
    1278             :    * if this is a level-1 element, this element's grandparent if this is
    1279             :    * a level-2 element, etc...
    1280             :    */
    1281             :   const Elem * top_parent () const;
    1282             : 
    1283             :   /**
    1284             :    * \returns The higher-dimensional Elem for which this Elem is a face.
    1285             :    *
    1286             :    * In some cases it is desirable to extract the boundary (or a subset thereof)
    1287             :    * of a D-dimensional mesh as a (D-1)-dimensional manifold.  In this case
    1288             :    * we may want to know the 'parent' element from which the manifold elements
    1289             :    * were extracted.  We can easily do that for the level-0 manifold elements
    1290             :    * by storing the D-dimensional parent.  This method provides access to that
    1291             :    * element.
    1292             :    *
    1293             :    * This method returns nullptr if this->dim() == LIBMESH_DIM; in
    1294             :    * such cases no data storage for an interior parent pointer has
    1295             :    * been allocated.
    1296             :    */
    1297             :   const Elem * interior_parent () const;
    1298             : 
    1299             :   Elem * interior_parent ();
    1300             : 
    1301             :   /**
    1302             :    * Sets the pointer to the element's interior_parent.
    1303             :    * Dangerous! Only use this if you know what you are doing!
    1304             :    */
    1305             :   void set_interior_parent (Elem * p);
    1306             : 
    1307             :   /**
    1308             :    * \returns The distance between nodes n1 and n2.
    1309             :    *
    1310             :    * Useful for computing the lengths of the sides of elements.
    1311             :    */
    1312             :   Real length (const unsigned int n1,
    1313             :                const unsigned int n2) const;
    1314             : 
    1315             :   /**
    1316             :    * \returns The number of adjacent vertices that uniquely define the
    1317             :    * location of the \f$ n^{th} \f$ second-order node, or 0 for linear
    1318             :    * elements.
    1319             :    *
    1320             :    * This method is useful when converting linear elements to quadratic
    1321             :    * elements.
    1322             :    *
    1323             :    * \note \p n has to be greater than or equal to \p this->n_vertices().
    1324             :    */
    1325             :   virtual unsigned int n_second_order_adjacent_vertices (const unsigned int n) const;
    1326             : 
    1327             :   /**
    1328             :    * \returns The element-local number of the \f$ v^{th} \f$ vertex
    1329             :    * that defines the \f$ n^{th} \f$ second-order node, or 0 for
    1330             :    * linear elements.
    1331             :    *
    1332             :    * \note The value is always less than \p this->n_vertices(), while
    1333             :    * \p n has to be greater than or equal to \p this->n_vertices().
    1334             :    */
    1335             :   virtual unsigned short int second_order_adjacent_vertex (const unsigned int n,
    1336             :                                                            const unsigned int v) const;
    1337             : 
    1338             :   /**
    1339             :    * \returns A pair (c,v), where
    1340             :    * c == child index, and
    1341             :    * v == element-local index of the \p \f$ n^{th} \f$
    1342             :    *      second-order node on the parent element.
    1343             :    * For linear elements, (0,0) is returned.
    1344             :    *
    1345             :    * \note The return values are always less than \p this->n_children()
    1346             :    * and \p this->child_ptr(c)->n_vertices().
    1347             :    *
    1348             :    * \note \p n has to be greater than or equal to \p this->n_vertices().
    1349             :    *
    1350             :    * \note On refined second-order elements, the return value will
    1351             :    * satisfy \p this->node_ptr(n) == this->child_ptr(c)->node_ptr(v).
    1352             :    */
    1353             :   virtual std::pair<unsigned short int, unsigned short int>
    1354             :   second_order_child_vertex (const unsigned int n) const;
    1355             : 
    1356             :   /**
    1357             :    * \returns The ElemType of the associated second-order element
    1358             :    * (which will be the same as the input if the input is already a
    1359             :    * second-order ElemType) or INVALID_ELEM for elements that cannot be
    1360             :    * converted into higher order equivalents.
    1361             :    *
    1362             :    * For example, when \p this is a \p TET4, then \p TET10 is returned.
    1363             :    *
    1364             :    * For some elements, there exist two second-order equivalents, e.g.
    1365             :    * for \p Quad4 there is \p Quad8 and \p Quad9.  When the optional
    1366             :    * \p full_ordered is \p true, then \p QUAD9 is returned.  When
    1367             :    * \p full_ordered is \p false, then \p QUAD8 is returned.
    1368             :    */
    1369             :   static ElemType second_order_equivalent_type (const ElemType et,
    1370             :                                                 const bool full_ordered=true);
    1371             : 
    1372             :   /**
    1373             :    * \returns The element type of the associated first-order element,
    1374             :    * or \p INVALID_ELEM for first-order or other elements that cannot be
    1375             :    * converted into lower order equivalents.
    1376             :    *
    1377             :    * For example, when \p this is a \p TET10, then \p TET4 is returned.
    1378             :    */
    1379             :   static ElemType first_order_equivalent_type (const ElemType et);
    1380             : 
    1381             :   /**
    1382             :    * \returns The ElemType of the associated "complete" order element
    1383             :    * (which will be the same as the input if the input is already a
    1384             :    * complete-order ElemType), or INVALID_ELEM for elements that cannot be
    1385             :    * converted into complete-order equivalents.
    1386             :    *
    1387             :    * The "complete" version of an element is an element which can
    1388             :    * represent the same geometry but which has nodes available to
    1389             :    * restore degrees of freedom on any vertex, edge, or face.
    1390             :    *
    1391             :    * For example, when \p this is a \p TET4, then \p TET14 is returned.
    1392             :    */
    1393             :   static ElemType complete_order_equivalent_type (const ElemType et);
    1394             : 
    1395             :   /**
    1396             :    * \returns The refinement level of the current element.
    1397             :    *
    1398             :    * If the element's parent is \p nullptr then by convention it is at
    1399             :    * level 0, otherwise it is simply at one level greater than its
    1400             :    * parent.
    1401             :    */
    1402             :   unsigned int level () const;
    1403             : 
    1404             :   /**
    1405             :    * \returns The value of the p refinement level of an active
    1406             :    * element, or the minimum value of the p refinement levels
    1407             :    * of an ancestor element's descendants.
    1408             :    */
    1409             :   unsigned int p_level () const;
    1410             : 
    1411             :   /**
    1412             :    * \returns \p true if the specified child is on the specified side.
    1413             :    */
    1414             :   virtual bool is_child_on_side(const unsigned int c,
    1415             :                                 const unsigned int s) const = 0;
    1416             : 
    1417             :   /**
    1418             :    * \returns The value of the mapping type for the element.
    1419             :    */
    1420             :   ElemMappingType mapping_type () const;
    1421             : 
    1422             :   /**
    1423             :    * Sets the value of the mapping type for the element.
    1424             :    */
    1425             :   void set_mapping_type (const ElemMappingType type);
    1426             : 
    1427             :   /**
    1428             :    * \returns The value of the mapping data for the element.
    1429             :    */
    1430             :   unsigned char mapping_data () const;
    1431             : 
    1432             :   /**
    1433             :    * Sets the value of the mapping data for the element.
    1434             :    */
    1435             :   void set_mapping_data (const unsigned char data);
    1436             : 
    1437             : 
    1438             : #ifdef LIBMESH_ENABLE_AMR
    1439             : 
    1440             :   /**
    1441             :    * Enumeration of possible element refinement states.
    1442             :    */
    1443             :   enum RefinementState { COARSEN = 0,
    1444             :                          DO_NOTHING,
    1445             :                          REFINE,
    1446             :                          JUST_REFINED,
    1447             :                          JUST_COARSENED,
    1448             :                          INACTIVE,
    1449             :                          COARSEN_INACTIVE,
    1450             :                          INVALID_REFINEMENTSTATE };
    1451             : 
    1452             :   /**
    1453             :    * \returns A constant pointer to the \f$ i^{th} \f$ child for this element.
    1454             :    * For internal use only - skips assertions about null pointers.
    1455             :    */
    1456             :   const Elem * raw_child_ptr (unsigned int i) const;
    1457             : 
    1458             :   /**
    1459             :    * \returns A constant pointer to the \f$ i^{th} \f$ child for this element.
    1460             :    * Do not call if this element has no children, i.e. is active.
    1461             :    */
    1462             :   const Elem * child_ptr (unsigned int i) const;
    1463             : 
    1464             :   /**
    1465             :    * \returns A non-constant pointer to the \f$ i^{th} \f$ child for this element.
    1466             :    * Do not call if this element has no children, i.e. is active.
    1467             :    */
    1468             :   Elem * child_ptr (unsigned int i);
    1469             : 
    1470             :   /**
    1471             :    * Nested classes for use iterating over all children of a parent
    1472             :    * element.
    1473             :    */
    1474             :   class ChildRefIter;
    1475             :   class ConstChildRefIter;
    1476             : 
    1477             :   /**
    1478             :    * Returns a range with all children of a parent element, usable in
    1479             :    * range-based for loops.  The exact type of the return value here
    1480             :    * may be subject to change in future libMesh releases, but the
    1481             :    * iterators will always dereference to produce a reference to a
    1482             :    * child element.
    1483             :    */
    1484             :   SimpleRange<ChildRefIter> child_ref_range();
    1485             : 
    1486             :   SimpleRange<ConstChildRefIter> child_ref_range() const;
    1487             : 
    1488             : private:
    1489             :   /**
    1490             :    * Sets the pointer to the \f$ i^{th} \f$ child for this element.
    1491             :    * Do not call if this element has no children, i.e. is active.
    1492             :    */
    1493             :   void set_child (unsigned int c, Elem * elem);
    1494             : 
    1495             : public:
    1496             :   /**
    1497             :    * \returns The child index which \p e corresponds to.
    1498             :    *
    1499             :    * I.e. if c = a->which_child_am_i(e); then a->child_ptr(c) will be
    1500             :    * e.
    1501             :    */
    1502             :   unsigned int which_child_am_i(const Elem * e) const;
    1503             : 
    1504             :   /**
    1505             :    * \returns \p true if the specified child is on the specified edge.
    1506             :    */
    1507             :   virtual bool is_child_on_edge(const unsigned int c,
    1508             :                                 const unsigned int e) const;
    1509             : 
    1510             :   /**
    1511             :    * Adds a child pointer to the array of children of this element.
    1512             :    * If this is the first child to be added, this method allocates
    1513             :    * memory in the parent's _children array, otherwise, it just sets
    1514             :    * the pointer.
    1515             :    */
    1516             :   void add_child (Elem * elem);
    1517             : 
    1518             :   /**
    1519             :    * Adds a new child pointer to the specified index in the array of
    1520             :    * children of this element.  If this is the first child to be added,
    1521             :    * this method allocates memory in the parent's _children array,
    1522             :    * otherwise, it just sets the pointer.
    1523             :    */
    1524             :   void add_child (Elem * elem, unsigned int c);
    1525             : 
    1526             :   /**
    1527             :    * Replaces the child pointer at the specified index in the child array.
    1528             :    */
    1529             :   void replace_child (Elem * elem, unsigned int c);
    1530             : 
    1531             :   /**
    1532             :    * Fills the vector \p family with the children of this element,
    1533             :    * recursively.  Calling this method on a twice-refined element
    1534             :    * will give you the element itself, its direct children, and their
    1535             :    * children, etc...  When the optional parameter \p reset is
    1536             :    * true, the vector will be cleared before the element and its
    1537             :    * descendants are added.
    1538             :    *
    1539             :    * The family tree only includes ancestor and active elements. To
    1540             :    * include subactive elements as well, use total_family_tree().
    1541             :    */
    1542             :   void family_tree (std::vector<const Elem *> & family,
    1543             :                     bool reset = true) const;
    1544             : 
    1545             :   /**
    1546             :    * Non-const version of function above; fills a vector of non-const pointers.
    1547             :    */
    1548             :   void family_tree (std::vector<Elem *> & family,
    1549             :                     bool reset = true);
    1550             : 
    1551             :   /**
    1552             :    * Same as the \p family_tree() member, but also adds any subactive
    1553             :    * descendants.
    1554             :    */
    1555             :   void total_family_tree (std::vector<const Elem *> & family,
    1556             :                           bool reset = true) const;
    1557             : 
    1558             :   /**
    1559             :    * Non-const version of function above; fills a vector of non-const pointers.
    1560             :    */
    1561             :   void total_family_tree (std::vector<Elem *> & family,
    1562             :                           bool reset = true);
    1563             : 
    1564             :   /**
    1565             :    * Same as the \p family_tree() member, but only adds the active
    1566             :    * children.  Can be thought of as removing all the inactive
    1567             :    * elements from the vector created by \p family_tree, but is
    1568             :    * implemented more efficiently.
    1569             :    */
    1570             :   void active_family_tree (std::vector<const Elem *> & active_family,
    1571             :                            bool reset = true) const;
    1572             : 
    1573             :   /**
    1574             :    * Non-const version of function above; fills a vector of non-const pointers.
    1575             :    */
    1576             :   void active_family_tree (std::vector<Elem *> & active_family,
    1577             :                            bool reset = true);
    1578             : 
    1579             :   /**
    1580             :    * Same as the \p family_tree() member, but only adds elements
    1581             :    * which are next to \p side.
    1582             :    */
    1583             :   void family_tree_by_side (std::vector<const Elem *> & family,
    1584             :                             unsigned int side,
    1585             :                             bool reset = true) const;
    1586             : 
    1587             :   /**
    1588             :    * Non-const version of function above; fills a vector of non-const pointers.
    1589             :    */
    1590             :   void family_tree_by_side (std::vector<Elem *> & family,
    1591             :                             unsigned int side,
    1592             :                             bool reset = true);
    1593             : 
    1594             :   /**
    1595             :    * Same as the \p total_family_tree() member, but only adds elements
    1596             :    * which are next to \p side.
    1597             :    */
    1598             :   void total_family_tree_by_side (std::vector<const Elem *> & family,
    1599             :                                   unsigned int side,
    1600             :                                   bool reset = true) const;
    1601             : 
    1602             :   /**
    1603             :    * Non-const version of function above; fills a vector of non-const pointers.
    1604             :    */
    1605             :   void total_family_tree_by_side (std::vector<Elem *> & family,
    1606             :                                   unsigned int side,
    1607             :                                   bool reset = true);
    1608             : 
    1609             :   /**
    1610             :    * Same as the \p active_family_tree() member, but only adds elements
    1611             :    * which are next to \p side.
    1612             :    */
    1613             :   void active_family_tree_by_side (std::vector<const Elem *> & family,
    1614             :                                    unsigned int side,
    1615             :                                    bool reset = true) const;
    1616             : 
    1617             :   /**
    1618             :    * Non-const version of function above; fills a vector of non-const pointers.
    1619             :    */
    1620             :   void active_family_tree_by_side (std::vector<Elem *> & family,
    1621             :                                    unsigned int side,
    1622             :                                    bool reset = true);
    1623             : 
    1624             :   /**
    1625             :    * Same as the \p family_tree() member, but only adds elements
    1626             :    * which are next to \p neighbor.
    1627             :    */
    1628             :   void family_tree_by_neighbor (std::vector<const Elem *> & family,
    1629             :                                 const Elem * neighbor,
    1630             :                                 bool reset = true) const;
    1631             : 
    1632             :   /**
    1633             :    * Non-const version of function above; fills a vector of non-const pointers.
    1634             :    */
    1635             :   void family_tree_by_neighbor (std::vector<Elem *> & family,
    1636             :                                 Elem * neighbor,
    1637             :                                 bool reset = true);
    1638             : 
    1639             :   /**
    1640             :    * Same as the \p family_tree_by_neighbor() member, but also adds
    1641             :    * any subactive descendants.
    1642             :    */
    1643             :   void total_family_tree_by_neighbor (std::vector<const Elem *> & family,
    1644             :                                       const Elem * neighbor,
    1645             :                                       bool reset = true) const;
    1646             : 
    1647             :   /**
    1648             :    * Non-const version of function above; fills a vector of non-const pointers.
    1649             :    */
    1650             :   void total_family_tree_by_neighbor (std::vector<Elem *> & family,
    1651             :                                       Elem * neighbor,
    1652             :                                       bool reset = true);
    1653             : 
    1654             :   /**
    1655             :    * Same as the \p family_tree() member, but only adds elements
    1656             :    * which are next to \p subneighbor.  Only applicable when
    1657             :    * \p this->has_neighbor(neighbor) and
    1658             :    * \p neighbor->is_ancestor(subneighbor)
    1659             :    */
    1660             :   void family_tree_by_subneighbor (std::vector<const Elem *> & family,
    1661             :                                    const Elem * neighbor,
    1662             :                                    const Elem * subneighbor,
    1663             :                                    bool reset = true) const;
    1664             : 
    1665             :   /**
    1666             :    * Non-const version of function above; fills a vector of non-const pointers.
    1667             :    */
    1668             :   void family_tree_by_subneighbor (std::vector<Elem *> & family,
    1669             :                                    Elem * neighbor,
    1670             :                                    Elem * subneighbor,
    1671             :                                    bool reset = true);
    1672             : 
    1673             :   /**
    1674             :    * Same as the \p family_tree_by_subneighbor() member, but also adds
    1675             :    * any subactive descendants.
    1676             :    */
    1677             :   void total_family_tree_by_subneighbor (std::vector<const Elem *> & family,
    1678             :                                          const Elem * neighbor,
    1679             :                                          const Elem * subneighbor,
    1680             :                                          bool reset = true) const;
    1681             : 
    1682             :   /**
    1683             :    * Non-const version of function above; fills a vector of non-const pointers.
    1684             :    */
    1685             :   void total_family_tree_by_subneighbor (std::vector<Elem *> & family,
    1686             :                                          Elem * neighbor,
    1687             :                                          Elem * subneighbor,
    1688             :                                          bool reset = true);
    1689             : 
    1690             :   /**
    1691             :    * Same as the \p active_family_tree() member, but only adds elements
    1692             :    * which are next to \p neighbor.
    1693             :    */
    1694             :   void active_family_tree_by_neighbor (std::vector<const Elem *> & family,
    1695             :                                        const Elem * neighbor,
    1696             :                                        bool reset = true) const;
    1697             : 
    1698             :   /**
    1699             :    * Non-const version of function above; fills a vector of non-const pointers.
    1700             :    */
    1701             :   void active_family_tree_by_neighbor (std::vector<Elem *> & family,
    1702             :                                        Elem * neighbor,
    1703             :                                        bool reset = true);
    1704             : 
    1705             :   /**
    1706             :    * Same as the \p active_family_tree_by_neighbor() member, but the
    1707             :    * \p neighbor here may be a topological (e.g. periodic boundary
    1708             :    * condition) neighbor, not just a local neighbor.
    1709             :    */
    1710             :   void active_family_tree_by_topological_neighbor (std::vector<const Elem *> & family,
    1711             :                                                    const Elem * neighbor,
    1712             :                                                    const MeshBase & mesh,
    1713             :                                                    const PointLocatorBase & point_locator,
    1714             :                                                    const PeriodicBoundaries * pb,
    1715             :                                                    bool reset = true) const;
    1716             : 
    1717             :   /**
    1718             :    * Non-const version of function above; fills a vector of non-const pointers.
    1719             :    */
    1720             :   void active_family_tree_by_topological_neighbor (std::vector<Elem *> & family,
    1721             :                                                    Elem * neighbor,
    1722             :                                                    const MeshBase & mesh,
    1723             :                                                    const PointLocatorBase & point_locator,
    1724             :                                                    const PeriodicBoundaries * pb,
    1725             :                                                    bool reset = true);
    1726             : 
    1727             :   /**
    1728             :    * \returns The value of the refinement flag for the element.
    1729             :    */
    1730             :   RefinementState refinement_flag () const;
    1731             : 
    1732             :   /**
    1733             :    * Sets the value of the refinement flag for the element.
    1734             :    */
    1735             :   void set_refinement_flag (const RefinementState rflag);
    1736             : 
    1737             :   /**
    1738             :    * \returns The value of the p-refinement flag for the element.
    1739             :    */
    1740             :   RefinementState p_refinement_flag () const;
    1741             : 
    1742             :   /**
    1743             :    * Sets the value of the p-refinement flag for the element.
    1744             :    */
    1745             :   void set_p_refinement_flag (const RefinementState pflag);
    1746             : 
    1747             :   /**
    1748             :    * \returns The maximum value of the p-refinement levels of
    1749             :    * an ancestor element's descendants.
    1750             :    */
    1751             :   unsigned int max_descendant_p_level () const;
    1752             : 
    1753             :   /**
    1754             :    * \returns The minimum p-refinement level of elements which are
    1755             :    * descended from this element, and which share a side with the
    1756             :    * active \p neighbor.
    1757             :    */
    1758             :   unsigned int min_p_level_by_neighbor (const Elem * neighbor,
    1759             :                                         unsigned int current_min) const;
    1760             : 
    1761             :   /**
    1762             :    * \returns The minimum new p-refinement level (i.e. after refinement
    1763             :    * and coarsening is done) of elements which are descended from this
    1764             :    * element and which share a side with the active \p neighbor.
    1765             :    */
    1766             :   unsigned int min_new_p_level_by_neighbor (const Elem * neighbor,
    1767             :                                             unsigned int current_min) const;
    1768             : 
    1769             :   /**
    1770             :    * Sets the value of the p-refinement level for the element.
    1771             :    *
    1772             :    * \note The maximum p-refinement level is currently 255.
    1773             :    */
    1774             :   void set_p_level (const unsigned int p);
    1775             : 
    1776             :   /**
    1777             :    * Sets the value of the p-refinement level for the element
    1778             :    * without altering the p-level of its ancestors
    1779             :    */
    1780             :   void hack_p_level (const unsigned int p);
    1781             : 
    1782             :   /**
    1783             :    * Sets the value of the p-refinement level for the element
    1784             :    * without altering the p-level of its ancestors; also sets the
    1785             :    * p_refinement_flag, simultaneously so that they can be safely
    1786             :    * checked for mutual consistency
    1787             :    */
    1788             :   void hack_p_level_and_refinement_flag (const unsigned int p,
    1789             :                                          RefinementState pflag);
    1790             : 
    1791             :   /**
    1792             :    * Refine the element.
    1793             :    */
    1794             :   virtual void refine (MeshRefinement & mesh_refinement);
    1795             : 
    1796             :   /**
    1797             :    * Coarsen the element.  This function is non-virtual since it is the same
    1798             :    * for all element types.
    1799             :    */
    1800             :   void coarsen ();
    1801             : 
    1802             :   /**
    1803             :    * Contract an active element, i.e. remove pointers to any
    1804             :    * subactive children.  This should only be called via
    1805             :    * MeshRefinement::contract, which will also remove subactive
    1806             :    * children from the mesh.
    1807             :    */
    1808             :   void contract ();
    1809             : 
    1810             : #endif
    1811             : 
    1812             : #ifndef NDEBUG
    1813             :   /**
    1814             :    * Checks for consistent neighbor links on this element.
    1815             :    */
    1816             :   void libmesh_assert_valid_neighbors() const;
    1817             : 
    1818             :   /**
    1819             :    * Checks for a valid id and pointers to nodes with valid ids on
    1820             :    * this element.
    1821             :    */
    1822             :   void libmesh_assert_valid_node_pointers() const;
    1823             : #endif // !NDEBUG
    1824             : 
    1825             :   /**
    1826             :    * \returns The local node index of the given point IF said node
    1827             :    * has a singular Jacobian for this element. If the given point
    1828             :    * is not a node or is a node and does not have a singular Jacobian,
    1829             :    * this will return invalid_uint.
    1830             :    *
    1831             :    * The intention is for this to be overridden in derived element
    1832             :    * classes that do have nodes that have singular Jacobians. When
    1833             :    * mapping failures are caught, we can check this to see if the
    1834             :    * failed physical point is actually a singular point and
    1835             :    * return the correct master point.
    1836             :    */
    1837           0 :   virtual unsigned int local_singular_node(const Point & /* p */, const Real /* tol */ = TOLERANCE*TOLERANCE) const
    1838           0 :   { return invalid_uint; }
    1839             : 
    1840             :   /**
    1841             :    * \returns true iff the node at the given index has a singular
    1842             :    * mapping; i.e. is the degree-4 node on a Pyramid.
    1843             :    */
    1844           0 :   virtual bool is_singular_node(unsigned int /* node_i */) const { return false; }
    1845             : 
    1846             :   /**
    1847             :    * \returns The local index of the center node on the side \p side.
    1848             :    *
    1849             :    * A center node is a node that is located at the centroid of the given side.
    1850             :    * If the given side does not have a center node, this will return invalid_uint.
    1851             :    */
    1852             :   virtual unsigned int center_node_on_side(const unsigned short side) const;
    1853             : 
    1854             : protected:
    1855             : 
    1856             :   /**
    1857             :    * The protected nested SideIter class is used to iterate over the
    1858             :    * sides of this Elem.  It is a specially-designed class since
    1859             :    * no sides are actually stored by the element.  This iterator-like
    1860             :    * class has to provide the following three operations
    1861             :    * 1) operator*
    1862             :    * 2) operator++
    1863             :    * 3) operator==
    1864             :    * The definition can be found at the end of this header file.
    1865             :    */
    1866             :   class SideIter;
    1867             : 
    1868             : public:
    1869             :   /**
    1870             :    * Useful iterator typedefs
    1871             :    */
    1872             :   typedef Predicates::multi_predicate Predicate;
    1873             : 
    1874             :   /**
    1875             :    * Data structure for iterating over sides.  Defined at the end of
    1876             :    * this header file.
    1877             :    */
    1878             :   struct side_iterator;
    1879             : 
    1880             :   /**
    1881             :    * Iterator accessor functions
    1882             :    */
    1883             :   side_iterator boundary_sides_begin();
    1884             :   side_iterator boundary_sides_end();
    1885             : 
    1886             : private:
    1887             :   /**
    1888             :    * Side iterator helper functions.  Used to replace the begin()
    1889             :    * and end() functions of the STL containers.
    1890             :    */
    1891             :   SideIter _first_side();
    1892             :   SideIter _last_side();
    1893             : 
    1894             : public:
    1895             : 
    1896             : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
    1897             : 
    1898             :   /**
    1899             :    * \returns \p true if the element is an infinite element,
    1900             :    * \p false otherwise.
    1901             :    */
    1902             :   virtual bool infinite () const = 0;
    1903             : 
    1904             :   /**
    1905             :    * \returns \p true if the specified (local) node number is a
    1906             :    * "mid-edge" node on an infinite element edge.
    1907             :    *
    1908             :    * This is false for all nodes on non-infinite elements, so we won't
    1909             :    * make it pure virtual, to simplify their code.
    1910             :    */
    1911           0 :   virtual bool is_mid_infinite_edge_node(const unsigned int /* n */) const
    1912           0 :   { libmesh_assert (!this->infinite()); return false; }
    1913             : 
    1914             :   /**
    1915             :    * \returns The origin for an infinite element.
    1916             :    *
    1917             :    * Currently, all infinite elements used in a mesh share the same
    1918             :    * origin.  Override this in infinite element classes.
    1919             :    */
    1920           0 :   virtual Point origin () const { libmesh_not_implemented(); return Point(); }
    1921             : 
    1922             : #else
    1923             : 
    1924             :   static constexpr bool infinite () { return false; }
    1925             : 
    1926             : #endif
    1927             : 
    1928             :   /**
    1929             :    * \returns An Elem of type \p type wrapped in a smart pointer.
    1930             :    */
    1931             :   static std::unique_ptr<Elem> build (const ElemType type,
    1932             :                                       Elem * p=nullptr);
    1933             : 
    1934             :   /**
    1935             :    * Calls the build() method above with a nullptr parent, and
    1936             :    * additionally sets the newly-created Elem's id. This can be useful
    1937             :    * when adding pre-numbered Elems to a Mesh via add_elem() calls.
    1938             :    */
    1939             :   static std::unique_ptr<Elem> build_with_id (const ElemType type,
    1940             :                                               dof_id_type id);
    1941             : 
    1942             :   /**
    1943             :    * \returns An Elem of the same type as \p this, wrapped in a smart
    1944             :    * pointer.
    1945             :    *
    1946             :    * This is not a complete clone() method (since e.g. it does not set
    1947             :    * node pointers; the standard use case reassigns node pointers from
    1948             :    * a different mesh), but it is necessary to use this instead of
    1949             :    * build() for runtime-polymorphic elements like Polygon subtypes
    1950             :    * whose "type" depends on more than their type(), and it is useful
    1951             :    * to use this for elements whose id, unique_id, extra integers,
    1952             :    * etc. should be preserved in the near-clone.
    1953             :    */
    1954             :   virtual std::unique_ptr<Elem> disconnected_clone () const;
    1955             : 
    1956             :   /**
    1957             :    * Returns the number of independent permutations of element nodes -
    1958             :    * e.g. a cube can be reoriented to put side 0 where side N is (for
    1959             :    * 0 <= N < 6) and then rotated in one of four ways, giving 24
    1960             :    * possible permutations.
    1961             :    *
    1962             :    * Permutations which change the mapping Jacobian of an element
    1963             :    * (i.e. flipping the element) are not allowed in this definition.
    1964             :    */
    1965             :   virtual unsigned int n_permutations() const = 0;
    1966             : 
    1967             :   /**
    1968             :    * Permutes the element (by swapping node and neighbor pointers)
    1969             :    * according to the specified index.
    1970             :    *
    1971             :    * This is useful for regression testing, by making it easy to make
    1972             :    * a structured mesh behave more like an arbitrarily unstructured
    1973             :    * mesh.
    1974             :    *
    1975             :    * This is so far *only* used for regression testing, so we do
    1976             :    * not currently provide a way to permute any boundary side/edge ids
    1977             :    * along with the element permutation.
    1978             :    */
    1979             :   virtual void permute(unsigned int perm_num) = 0;
    1980             : 
    1981             :   /**
    1982             :    * Flips the element (by swapping node and neighbor pointers) to
    1983             :    * have a mapping Jacobian of opposite sign.
    1984             :    *
    1985             :    * This is useful for automatically fixing up elements that have
    1986             :    * been newly created (e.g. from extrusions) with a negative
    1987             :    * Jacobian.
    1988             :    *
    1989             :    * If \p boundary_info is not null, swap boundary side/edge ids
    1990             :    * consistently.
    1991             :    */
    1992             :   virtual void flip(BoundaryInfo * boundary_info) = 0;
    1993             : 
    1994             :   /**
    1995             :    * \returns Whether the element is flipped compared to standard
    1996             :    * libMesh (e.g. clockwise for 2D elements) node orientations.
    1997             :    *
    1998             :    * Always returns \p false if a 2D element is not in the XY plane or
    1999             :    * a 1D element is not on the X axis; user code designed to work for
    2000             :    * embedded manifolds should handle any consistent orientation, and
    2001             :    * determining whether an orientation is consistent is not a local
    2002             :    * operation.
    2003             :    */
    2004             :    virtual bool is_flipped() const = 0;
    2005             : 
    2006             :   /**
    2007             :    * Flips the element (by swapping node and neighbor pointers) to
    2008             :    * have a mapping Jacobian of opposite sign, iff we find a negative
    2009             :    * orientation.  This only fixes flipped elements; for tangled
    2010             :    * elements the only fixes possible are non-local.
    2011             :    */
    2012             :   void orient(BoundaryInfo * boundary_info);
    2013             : 
    2014             : #ifdef LIBMESH_ENABLE_AMR
    2015             : 
    2016             :   /**
    2017             :    * \returns The local node id on the parent which corresponds to node
    2018             :    * \p n of child \p c, or \p invalid_uint if no such parent
    2019             :    * node exists.
    2020             :    */
    2021             :   virtual unsigned int as_parent_node (unsigned int c,
    2022             :                                        unsigned int n) const;
    2023             : 
    2024             :   /**
    2025             :    * \returns All the pairs of nodes (indexed by local node id) which
    2026             :    * should bracket node \p n of child \p c.
    2027             :    */
    2028             :   virtual
    2029             :   const std::vector<std::pair<unsigned char, unsigned char>> &
    2030             :   parent_bracketing_nodes(unsigned int c,
    2031             :                           unsigned int n) const;
    2032             : 
    2033             :   /**
    2034             :    * \returns All the pairs of nodes (indexed by global node id) which
    2035             :    * should bracket node \p n of child \p c.
    2036             :    */
    2037             :   virtual
    2038             :   const std::vector<std::pair<dof_id_type, dof_id_type>>
    2039             :   bracketing_nodes(unsigned int c,
    2040             :                    unsigned int n) const;
    2041             : 
    2042             : 
    2043             :   /**
    2044             :    * \returns The embedding matrix entry for the requested child.
    2045             :    */
    2046             :   virtual Real embedding_matrix (const unsigned int child_num,
    2047             :                                  const unsigned int child_node_num,
    2048             :                                  const unsigned int parent_node_num) const = 0;
    2049             : 
    2050             :   /**
    2051             :    * \returns A "version number" that identifies which embedding
    2052             :    * matrix is in use.
    2053             :    *
    2054             :    * Some element types may use a different embedding matrix depending
    2055             :    * on their geometric characteristics.
    2056             :    */
    2057           0 :   virtual unsigned int embedding_matrix_version () const { return 0; }
    2058             : 
    2059             : #endif // LIBMESH_ENABLE_AMR
    2060             : 
    2061             : 
    2062             : protected:
    2063             : 
    2064             :   /**
    2065             :    * Default tolerance to use in has_affine_map().
    2066             :    */
    2067             :   static constexpr Real affine_tol = TOLERANCE*TOLERANCE;
    2068             : 
    2069             :   /**
    2070             :    * \returns A hash key computed from a single node id.
    2071             :    */
    2072             :   static dof_id_type compute_key (dof_id_type n0);
    2073             : 
    2074             :   /**
    2075             :    * \returns A hash key computed from two node ids.
    2076             :    */
    2077             :   static dof_id_type compute_key (dof_id_type n0,
    2078             :                                   dof_id_type n1);
    2079             : 
    2080             :   /**
    2081             :    * \returns A hash key computed from three node ids.
    2082             :    */
    2083             :   static dof_id_type compute_key (dof_id_type n0,
    2084             :                                   dof_id_type n1,
    2085             :                                   dof_id_type n2);
    2086             : 
    2087             :   /**
    2088             :    * \returns A hash key computed from four node ids.
    2089             :    */
    2090             :   static dof_id_type compute_key (dof_id_type n0,
    2091             :                                   dof_id_type n1,
    2092             :                                   dof_id_type n2,
    2093             :                                   dof_id_type n3);
    2094             : 
    2095             :   /**
    2096             :    * Swaps two node_ptrs
    2097             :    */
    2098    27628562 :   void swap2nodes(unsigned int n1, unsigned int n2)
    2099             :   {
    2100     3399728 :     Node * temp = this->node_ptr(n1);
    2101    29328426 :     this->set_node(n1, this->node_ptr(n2));
    2102    27628562 :     this->set_node(n2, temp);
    2103    27628562 :   }
    2104             : 
    2105             :   /**
    2106             :    * Swaps two neighbor_ptrs
    2107             :    */
    2108     6631539 :   void swap2neighbors(unsigned int n1, unsigned int n2)
    2109             :   {
    2110      646486 :     Elem * temp = this->neighbor_ptr(n1);
    2111      444246 :     this->set_neighbor(n1, this->neighbor_ptr(n2));
    2112      444246 :     this->set_neighbor(n2, temp);
    2113     6705663 :   }
    2114             : 
    2115             :   /**
    2116             :    * Swaps two sides in \p boundary_info, if it is non-null.
    2117             :    */
    2118             :   void swap2boundarysides(unsigned short s1, unsigned short s2,
    2119             :                           BoundaryInfo * boundary_info) const;
    2120             : 
    2121             :   /**
    2122             :    * Swaps two edges in \p boundary_info, if it is non-null.
    2123             :    */
    2124             :   void swap2boundaryedges(unsigned short e1, unsigned short e2,
    2125             :                           BoundaryInfo * boundary_info) const;
    2126             : 
    2127             :   /**
    2128             :    * Swaps three node_ptrs, "rotating" them.
    2129             :    */
    2130    11044046 :   void swap3nodes(unsigned int n1, unsigned int n2, unsigned int n3)
    2131             :   {
    2132    11741906 :     swap2nodes(n1, n2);
    2133    11741906 :     swap2nodes(n2, n3);
    2134    11044046 :   }
    2135             : 
    2136             :   /**
    2137             :    * Swaps three neighbor_ptrs, "rotating" them.
    2138             :    */
    2139     3009388 :   void swap3neighbors(unsigned int n1, unsigned int n2,
    2140             :                       unsigned int n3)
    2141             :   {
    2142     3009388 :     swap2neighbors(n1, n2);
    2143     3009388 :     swap2neighbors(n2, n3);
    2144     3009388 :   }
    2145             : 
    2146             :   /**
    2147             :    * Swaps four node_ptrs, "rotating" them.
    2148             :    */
    2149     3036988 :   void swap4nodes(unsigned int n1, unsigned int n2, unsigned int n3,
    2150             :                   unsigned int n4)
    2151             :   {
    2152     2806960 :     swap3nodes(n1, n2, n3);
    2153     3036988 :     swap2nodes(n3, n4);
    2154     3036988 :   }
    2155             : 
    2156             :   /**
    2157             :    * Swaps four neighbor_ptrs, "rotating" them.
    2158             :    */
    2159      590067 :   void swap4neighbors(unsigned int n1, unsigned int n2,
    2160             :                       unsigned int n3, unsigned int n4)
    2161             :   {
    2162      590067 :     swap3neighbors(n1, n2, n3);
    2163      590067 :     swap2neighbors(n3, n4);
    2164      590067 :   }
    2165             : 
    2166             : 
    2167             :   /**
    2168             :    * An implementation for simple (all sides equal) elements
    2169             :    */
    2170             :   template <typename Sideclass, typename Subclass>
    2171             :   std::unique_ptr<Elem>
    2172             :   simple_build_side_ptr(const unsigned int i);
    2173             : 
    2174             :   /**
    2175             :    * An implementation for simple (all sides equal) elements
    2176             :    */
    2177             :   template <typename Subclass>
    2178             :   void simple_build_side_ptr(std::unique_ptr<Elem> & side,
    2179             :                              const unsigned int i,
    2180             :                              ElemType sidetype);
    2181             : 
    2182             :   /**
    2183             :    * An implementation for simple (all sides equal) elements
    2184             :    */
    2185             :   template <typename Subclass, typename Mapclass>
    2186             :   void simple_side_ptr(std::unique_ptr<Elem> & side,
    2187             :                        const unsigned int i,
    2188             :                        ElemType sidetype);
    2189             : 
    2190             :   /**
    2191             :    * An implementation for simple (all edges equal) elements
    2192             :    */
    2193             :   template <typename Edgeclass, typename Subclass>
    2194             :   std::unique_ptr<Elem>
    2195             :   simple_build_edge_ptr(const unsigned int i);
    2196             : 
    2197             :   /**
    2198             :    * An implementation for simple (all edges equal) elements
    2199             :    */
    2200             :   template <typename Subclass>
    2201             :   void simple_build_edge_ptr(std::unique_ptr<Elem> & edge,
    2202             :                              const unsigned int i,
    2203             :                              ElemType edgetype);
    2204             : 
    2205             : 
    2206             : #ifdef LIBMESH_ENABLE_AMR
    2207             : 
    2208             :   /**
    2209             :    * Elem subclasses which don't do their own bracketing node
    2210             :    * calculations will need to supply a static cache, since the
    2211             :    * default calculation is slow.
    2212             :    */
    2213             :   virtual
    2214             :   std::vector<std::vector<std::vector<std::vector<std::pair<unsigned char, unsigned char>>>>> &
    2215           0 :   _get_bracketing_node_cache() const
    2216             :   {
    2217           0 :     static std::vector<std::vector<std::vector<std::vector<std::pair<unsigned char, unsigned char>>>>> c;
    2218           0 :     libmesh_error();
    2219             :     return c;
    2220             :   }
    2221             : 
    2222             :   /**
    2223             :    * Elem subclasses which don't do their own child-to-parent node
    2224             :    * calculations will need to supply a static cache, since the
    2225             :    * default calculation is slow.
    2226             :    */
    2227             :   virtual
    2228             :   std::vector<std::vector<std::vector<signed char>>> &
    2229           0 :   _get_parent_indices_cache() const
    2230             :   {
    2231           0 :     static std::vector<std::vector<std::vector<signed char>>> c;
    2232           0 :     libmesh_error();
    2233             :     return c;
    2234             :   }
    2235             : 
    2236             : #endif // LIBMESH_ENABLE_AMR
    2237             : 
    2238             : public:
    2239             : 
    2240             :   /**
    2241             :    * Replaces this element with \p nullptr for all of its neighbors.
    2242             :    * This is useful when deleting an element.
    2243             :    */
    2244             :   void nullify_neighbors ();
    2245             : 
    2246             : protected:
    2247             : 
    2248             :   /**
    2249             :    * Pointers to the nodes we are connected to.
    2250             :    */
    2251             :   Node ** _nodes;
    2252             : 
    2253             :   /**
    2254             :    * Pointers to this element's parent and neighbors, and for
    2255             :    * lower-dimensional elements' interior_parent.
    2256             :    */
    2257             :   Elem ** _elemlinks;
    2258             : 
    2259             : #ifdef LIBMESH_ENABLE_AMR
    2260             :   /**
    2261             :    * unique_ptr to array of this element's children.
    2262             :    *
    2263             :    * A Mesh ultimately owns the child Elems so we are not responsible
    2264             :    * for deleting them, but we are responsible for cleaning up the
    2265             :    * array allocated to hold those Elems, hence the unique_ptr.
    2266             :    */
    2267             :   std::unique_ptr<Elem *[]> _children;
    2268             : #endif
    2269             : 
    2270             :   /**
    2271             :    * The subdomain to which this element belongs.
    2272             :    */
    2273             :   subdomain_id_type _sbd_id;
    2274             : 
    2275             : #ifdef LIBMESH_ENABLE_AMR
    2276             :   /**
    2277             :    * h refinement flag. This is stored as an unsigned char
    2278             :    * to save space.
    2279             :    */
    2280             :   unsigned char _rflag;
    2281             : 
    2282             :   /**
    2283             :    * p refinement flag. This is stored as an unsigned char
    2284             :    * to save space.
    2285             :    */
    2286             :   unsigned char _pflag;
    2287             : 
    2288             :   /**
    2289             :    * p refinement level - the difference between the
    2290             :    * polynomial degree on this element and the minimum
    2291             :    * polynomial degree on the mesh.
    2292             :    * This is stored as an unsigned char to save space.
    2293             :    * In theory, these last four bytes might have
    2294             :    * been padding anyway.
    2295             :    */
    2296             :   unsigned char _p_level;
    2297             : #endif
    2298             : 
    2299             :   /**
    2300             :    * Mapping function type; currently either 0 (LAGRANGE) or 1
    2301             :    * (RATIONAL_BERNSTEIN).
    2302             :    */
    2303             :   unsigned char _map_type;
    2304             : 
    2305             :   /**
    2306             :    * Mapping function data; currently used when needed to store the
    2307             :    * RATIONAL_BERNSTEIN nodal weight data index.
    2308             :    */
    2309             :   unsigned char _map_data;
    2310             : };
    2311             : 
    2312             : 
    2313             : 
    2314             : // ------------------------------------------------------------
    2315             : // Elem helper classes
    2316             : //
    2317             : class
    2318             : Elem::NodeRefIter : public PointerToPointerIter<Node>
    2319             : {
    2320             : public:
    2321    22703510 :   NodeRefIter (Node * const * nodepp) : PointerToPointerIter<Node>(nodepp) {}
    2322             : };
    2323             : 
    2324             : 
    2325             : class
    2326             : Elem::ConstNodeRefIter : public PointerToPointerIter<const Node>
    2327             : {
    2328             : public:
    2329    17734716 :   ConstNodeRefIter (const Node * const * nodepp) : PointerToPointerIter<const Node>(nodepp) {}
    2330             : };
    2331             : 
    2332             : 
    2333             : #ifdef LIBMESH_ENABLE_AMR
    2334             : class
    2335             : Elem::ChildRefIter : public PointerToPointerIter<Elem>
    2336             : {
    2337             : public:
    2338     6905820 :   ChildRefIter (Elem * const * childpp) : PointerToPointerIter<Elem>(childpp) {}
    2339             : };
    2340             : 
    2341             : 
    2342             : class
    2343             : Elem::ConstChildRefIter : public PointerToPointerIter<const Elem>
    2344             : {
    2345             : public:
    2346     2756128 :   ConstChildRefIter (const Elem * const * childpp) : PointerToPointerIter<const Elem>(childpp) {}
    2347             : };
    2348             : 
    2349             : 
    2350             : 
    2351             : inline
    2352    64749691 : SimpleRange<Elem::ChildRefIter> Elem::child_ref_range()
    2353             : {
    2354     3452910 :   libmesh_assert(_children);
    2355    67606791 :   return {_children.get(), _children.get() + this->n_children()};
    2356             : }
    2357             : 
    2358             : 
    2359             : inline
    2360     7501563 : SimpleRange<Elem::ConstChildRefIter> Elem::child_ref_range() const
    2361             : {
    2362     1378064 :   libmesh_assert(_children);
    2363     7787215 :   return {_children.get(), _children.get() + this->n_children()};
    2364             : }
    2365             : #endif // LIBMESH_ENABLE_AMR
    2366             : 
    2367             : 
    2368             : 
    2369             : 
    2370             : // ------------------------------------------------------------
    2371             : // global Elem functions
    2372             : 
    2373             : inline
    2374           0 : std::ostream & operator << (std::ostream & os, const Elem & e)
    2375             : {
    2376           0 :   e.print_info(os);
    2377           0 :   return os;
    2378             : }
    2379             : 
    2380             : 
    2381             : // ------------------------------------------------------------
    2382             : // Elem class member functions
    2383             : inline
    2384   502425391 : Elem::Elem(const unsigned int nn,
    2385             :            const unsigned int ns,
    2386             :            Elem * p,
    2387             :            Elem ** elemlinkdata,
    2388   502425391 :            Node ** nodelinkdata) :
    2389   125787380 :   _nodes(nodelinkdata),
    2390   125787380 :   _elemlinks(elemlinkdata),
    2391   125787380 :   _sbd_id(0),
    2392             : #ifdef LIBMESH_ENABLE_AMR
    2393   125787380 :   _rflag(Elem::DO_NOTHING),
    2394   125787380 :   _pflag(Elem::DO_NOTHING),
    2395   125787380 :   _p_level(0),
    2396             : #endif
    2397   126384892 :   _map_type(p ? p->mapping_type() : 0),
    2398   531512503 :   _map_data(p ? p->mapping_data() : 0)
    2399             : {
    2400   502425391 :   this->processor_id() = DofObject::invalid_processor_id;
    2401             : 
    2402             :   // If this ever legitimately fails we need to increase max_n_nodes
    2403    36724171 :   libmesh_assert_less_equal(nn, max_n_nodes);
    2404             : 
    2405             :   // We currently only support refinement of elements into child
    2406             :   // elements of the same type.  We can't test elem->type() here,
    2407             :   // because that's virtual and we're still in the base class
    2408             :   // constructor, but we can at least usually verify constency with
    2409             :   // the arguments we were handed.
    2410             : #ifndef NDEBUG
    2411    36724171 :   if (p && !p->runtime_topology())
    2412             :     {
    2413      294636 :       libmesh_assert_equal_to(nn, p->n_nodes());
    2414      294636 :       libmesh_assert_equal_to(ns, p->n_sides());
    2415             :     }
    2416             : #endif
    2417             : 
    2418             :   // Initialize the nodes data structure if we're given a pointer to
    2419             :   // memory for it.
    2420   502425391 :   if (_nodes)
    2421             :     {
    2422  2312052971 :       for (unsigned int n=0; n<nn; n++)
    2423  1809720506 :         _nodes[n] = nullptr;
    2424             :     }
    2425             : 
    2426             :   // Initialize the neighbors/parent data structure
    2427             :   // _elemlinks = new Elem *[ns+1];
    2428             : 
    2429             :   // Initialize the elements data structure if we're given a pointer
    2430             :   // to memory for it.  If we *weren't* given memory for it, e.g.
    2431             :   // because a subclass like an arbitrary Polygon needs to
    2432             :   // heap-allocate this memory, then that subclass will have to handle
    2433             :   // this initialization too.
    2434   502425391 :   if (_elemlinks)
    2435             :     {
    2436   502349063 :       _elemlinks[0] = p;
    2437             : 
    2438  2114024161 :       for (unsigned int n=1; n<ns+1; n++)
    2439  1611675098 :         _elemlinks[n] = nullptr;
    2440             : 
    2441             :       // Optionally initialize data from the parent
    2442   502349063 :       if (this->parent())
    2443             :         {
    2444    29087112 :           this->subdomain_id() = this->parent()->subdomain_id();
    2445    29087112 :           this->processor_id() = this->parent()->processor_id();
    2446    29087112 :           _map_type = this->parent()->_map_type;
    2447    29087112 :           _map_data = this->parent()->_map_data;
    2448             : 
    2449             : #ifdef LIBMESH_ENABLE_AMR
    2450    29389988 :           this->set_p_level(this->parent()->p_level());
    2451             : #endif
    2452             :         }
    2453             :     }
    2454   502425391 : }
    2455             : 
    2456             : 
    2457             : 
    2458             : inline
    2459  3468286049 : const Point & Elem::point (const unsigned int i) const
    2460             : {
    2461  3468286049 :   libmesh_assert_less (i, this->n_nodes());
    2462  3468286049 :   libmesh_assert(_nodes[i]);
    2463  3468286049 :   libmesh_assert_not_equal_to (_nodes[i]->id(), Node::invalid_id);
    2464             : 
    2465 51279299647 :   return *_nodes[i];
    2466             : }
    2467             : 
    2468             : 
    2469             : 
    2470             : inline
    2471     2871960 : Point & Elem::point (const unsigned int i)
    2472             : {
    2473     2871960 :   libmesh_assert_less (i, this->n_nodes());
    2474             : 
    2475   175028778 :   return *_nodes[i];
    2476             : }
    2477             : 
    2478             : 
    2479             : 
    2480             : inline
    2481    91367604 : dof_id_type Elem::node_id (const unsigned int i) const
    2482             : {
    2483    91367604 :   libmesh_assert_less (i, this->n_nodes());
    2484    91367604 :   libmesh_assert(_nodes[i]);
    2485    91367604 :   libmesh_assert_not_equal_to (_nodes[i]->id(), Node::invalid_id);
    2486             : 
    2487   997171422 :   return _nodes[i]->id();
    2488             : }
    2489             : 
    2490             : 
    2491             : 
    2492             : inline
    2493     1218369 : unsigned int Elem::local_node (const dof_id_type i) const
    2494             : {
    2495     5583811 :   for (auto n : make_range(this->n_nodes()))
    2496     5583811 :     if (this->node_id(n) == i)
    2497        2328 :       return n;
    2498             : 
    2499           0 :   return libMesh::invalid_uint;
    2500             : }
    2501             : 
    2502             : 
    2503             : 
    2504             : inline
    2505    25103292 : const Node * const * Elem::get_nodes () const
    2506             : {
    2507   137523291 :   return _nodes;
    2508             : }
    2509             : 
    2510             : 
    2511             : 
    2512             : inline
    2513   115236295 : const Node * Elem::node_ptr (const unsigned int i) const
    2514             : {
    2515   115236295 :   libmesh_assert_less (i, this->n_nodes());
    2516   115236295 :   libmesh_assert(_nodes[i]);
    2517             : 
    2518 25104298082 :   return _nodes[i];
    2519             : }
    2520             : 
    2521             : 
    2522             : 
    2523             : inline
    2524   129308113 : Node * Elem::node_ptr (const unsigned int i)
    2525             : {
    2526   129308113 :   libmesh_assert_less (i, this->n_nodes());
    2527   129308113 :   libmesh_assert(_nodes[i]);
    2528             : 
    2529  1481409545 :   return _nodes[i];
    2530             : }
    2531             : 
    2532             : 
    2533             : 
    2534             : inline
    2535    19913078 : const Node & Elem::node_ref (const unsigned int i) const
    2536             : {
    2537   185178413 :   return *this->node_ptr(i);
    2538             : }
    2539             : 
    2540             : 
    2541             : 
    2542             : inline
    2543    25689948 : Node & Elem::node_ref (const unsigned int i)
    2544             : {
    2545    53902032 :   return *this->node_ptr(i);
    2546             : }
    2547             : 
    2548             : 
    2549             : 
    2550             : inline
    2551      859366 : unsigned int Elem::get_node_index (const Node * node_ptr) const
    2552             : {
    2553     3041023 :   for (auto n : make_range(this->n_nodes()))
    2554     3026857 :     if (this->_nodes[n] == node_ptr)
    2555      161381 :       return n;
    2556             : 
    2557           0 :   return libMesh::invalid_uint;
    2558             : }
    2559             : 
    2560             : 
    2561             : 
    2562             : #ifdef LIBMESH_ENABLE_DEPRECATED
    2563             : inline
    2564           0 : Node * & Elem::set_node (const unsigned int i)
    2565             : {
    2566           0 :   libmesh_assert_less (i, this->n_nodes());
    2567             : 
    2568             :   libmesh_deprecated();
    2569             : 
    2570           0 :   return _nodes[i];
    2571             : }
    2572             : #endif // LIBMESH_ENABLE_DEPRECATED
    2573             : 
    2574             : 
    2575             : 
    2576             : inline
    2577  3098720587 : void Elem::set_node (const unsigned int i,
    2578             :                      Node * node)
    2579             : {
    2580      177702 :   libmesh_assert_less (i, this->n_nodes());
    2581             : 
    2582  3386910102 :   _nodes[i] = node;
    2583  2734802887 : }
    2584             : 
    2585             : 
    2586             : 
    2587             : inline
    2588    59997810 : subdomain_id_type Elem::subdomain_id () const
    2589             : {
    2590   553814503 :   return _sbd_id;
    2591             : }
    2592             : 
    2593             : 
    2594             : 
    2595             : inline
    2596    44433912 : subdomain_id_type & Elem::subdomain_id ()
    2597             : {
    2598   120265696 :   return _sbd_id;
    2599             : }
    2600             : 
    2601             : 
    2602             : 
    2603             : inline
    2604           0 : bool Elem::operator != (const Elem & rhs) const
    2605             : {
    2606       25041 :   return !(*this == rhs);
    2607             : }
    2608             : 
    2609             : 
    2610             : 
    2611             : inline
    2612   103750892 : const Elem * Elem::neighbor_ptr (unsigned int i) const
    2613             : {
    2614   103750892 :   libmesh_assert_less (i, this->n_neighbors());
    2615             : 
    2616   494641804 :   return _elemlinks[i+1];
    2617             : }
    2618             : 
    2619             : 
    2620             : 
    2621             : inline
    2622    18270223 : Elem * Elem::neighbor_ptr (unsigned int i)
    2623             : {
    2624    18270223 :   libmesh_assert_less (i, this->n_neighbors());
    2625             : 
    2626  1451915777 :   return _elemlinks[i+1];
    2627             : }
    2628             : 
    2629             : 
    2630             : 
    2631             : inline
    2632     9560978 : void Elem::set_neighbor (const unsigned int i, Elem * n)
    2633             : {
    2634     9560978 :   libmesh_assert_less (i, this->n_neighbors());
    2635             : 
    2636   577605354 :   _elemlinks[i+1] = n;
    2637   692835316 : }
    2638             : 
    2639             : 
    2640             : 
    2641             : inline
    2642    94496967 : bool Elem::has_neighbor (const Elem * elem) const
    2643             : {
    2644   283125840 :   for (auto n : this->neighbor_ptr_range())
    2645   271512466 :     if (n == elem)
    2646    83240130 :       return true;
    2647             : 
    2648    11256837 :   return false;
    2649             : }
    2650             : 
    2651             : 
    2652             : 
    2653             : inline
    2654             : Elem * Elem::child_neighbor (Elem * elem)
    2655             : {
    2656             :   for (auto n : elem->neighbor_ptr_range())
    2657             :     if (n && n->parent() == this)
    2658             :       return n;
    2659             : 
    2660             :   return nullptr;
    2661             : }
    2662             : 
    2663             : 
    2664             : 
    2665             : inline
    2666             : const Elem * Elem::child_neighbor (const Elem * elem) const
    2667             : {
    2668             :   for (auto n : elem->neighbor_ptr_range())
    2669             :     if (n && n->parent() == this)
    2670             :       return n;
    2671             : 
    2672             :   return nullptr;
    2673             : }
    2674             : 
    2675             : 
    2676             : 
    2677             : inline
    2678             : SimpleRange<Elem::NodeRefIter>
    2679   207753482 : Elem::node_ref_range()
    2680             : {
    2681   217426196 :   return {_nodes, _nodes+this->n_nodes()};
    2682             : }
    2683             : 
    2684             : 
    2685             : 
    2686             : inline
    2687             : SimpleRange<Elem::ConstNodeRefIter>
    2688   268171784 : Elem::node_ref_range() const
    2689             : {
    2690   285758276 :   return {_nodes, _nodes+this->n_nodes()};
    2691             : }
    2692             : 
    2693             : 
    2694             : 
    2695             : inline
    2696             : IntRange<unsigned short>
    2697    40562575 : Elem::node_index_range() const
    2698             : {
    2699   584535811 :   return {0, cast_int<unsigned short>(this->n_nodes())};
    2700             : }
    2701             : 
    2702             : 
    2703             : 
    2704             : inline
    2705             : IntRange<unsigned short>
    2706      384081 : Elem::edge_index_range() const
    2707             : {
    2708    48282185 :   return {0, cast_int<unsigned short>(this->n_edges())};
    2709             : }
    2710             : 
    2711             : 
    2712             : 
    2713             : inline
    2714             : IntRange<unsigned short>
    2715      293042 : Elem::face_index_range() const
    2716             : {
    2717     3888420 :   return {0, cast_int<unsigned short>(this->n_faces())};
    2718             : }
    2719             : 
    2720             : 
    2721             : 
    2722             : inline
    2723             : IntRange<unsigned short>
    2724     8727035 : Elem::side_index_range() const
    2725             : {
    2726   379597019 :   return {0, cast_int<unsigned short>(this->n_sides())};
    2727             : }
    2728             : 
    2729             : 
    2730             : 
    2731             : 
    2732             : inline
    2733             : std::unique_ptr<const Elem> Elem::side_ptr (unsigned int i) const
    2734             : {
    2735             :   // Call the non-const version of this function, return the result as
    2736             :   // a std::unique_ptr<const Elem>.
    2737             :   Elem * me = const_cast<Elem *>(this);
    2738             :   return me->side_ptr(i);
    2739             : }
    2740             : 
    2741             : 
    2742             : 
    2743             : inline
    2744             : void
    2745       15904 : Elem::side_ptr (std::unique_ptr<const Elem> & elem,
    2746             :                 const unsigned int i) const
    2747             : {
    2748             :   // Hand off to the non-const version of this function
    2749         448 :   Elem * me = const_cast<Elem *>(this);
    2750         896 :   std::unique_ptr<Elem> e {const_cast<Elem *>(elem.release())};
    2751       15904 :   me->side_ptr(e, i);
    2752       15456 :   elem = std::move(e);
    2753       15904 : }
    2754             : 
    2755             : 
    2756             : 
    2757             : inline
    2758             : std::unique_ptr<const Elem>
    2759   108202401 : Elem::build_side_ptr (const unsigned int i) const
    2760             : {
    2761             :   // Call the non-const version of this function, return the result as
    2762             :   // a std::unique_ptr<const Elem>.
    2763     2909904 :   Elem * me = const_cast<Elem *>(this);
    2764   118463204 :   return me->build_side_ptr(i);
    2765             : }
    2766             : 
    2767             : 
    2768             : 
    2769             : inline
    2770             : void
    2771    42114295 : Elem::build_side_ptr (std::unique_ptr<const Elem> & elem,
    2772             :                       const unsigned int i) const
    2773             : {
    2774             :   // Hand off to the non-const version of this function
    2775      755780 :   Elem * me = const_cast<Elem *>(this);
    2776     1511560 :   std::unique_ptr<Elem> e {const_cast<Elem *>(elem.release())};
    2777    42114295 :   me->build_side_ptr(e, i);
    2778    40416439 :   elem = std::move(e);
    2779    42114295 : }
    2780             : 
    2781             : 
    2782             : 
    2783             : template <typename Sideclass, typename Subclass>
    2784             : inline
    2785             : std::unique_ptr<Elem>
    2786   159791227 : Elem::simple_build_side_ptr (const unsigned int i)
    2787             : {
    2788    34925368 :   libmesh_assert_less (i, this->n_sides());
    2789             : 
    2790   159791227 :   std::unique_ptr<Elem> face = std::make_unique<Sideclass>();
    2791  1218689719 :   for (auto n : face->node_index_range())
    2792  1058898492 :     face->set_node(n, this->node_ptr(Subclass::side_nodes_map[i][n]));
    2793             : 
    2794   159791227 :   face->set_interior_parent(this);
    2795   148675967 :   face->inherit_data_from(*this);
    2796             : 
    2797   159791227 :   return face;
    2798           0 : }
    2799             : 
    2800             : 
    2801             : 
    2802             : template <typename Subclass>
    2803             : inline
    2804             : void
    2805    41325596 : Elem::simple_build_side_ptr (std::unique_ptr<Elem> & side,
    2806             :                              const unsigned int i,
    2807             :                              ElemType sidetype)
    2808             : {
    2809     1677553 :   libmesh_assert_less (i, this->n_sides());
    2810             : 
    2811    41325596 :   if (!side.get() || side->type() != sidetype)
    2812             :     {
    2813      172221 :       Subclass & real_me = cast_ref<Subclass&>(*this);
    2814     1540155 :       side = real_me.Subclass::build_side_ptr(i);
    2815             :     }
    2816             :   else
    2817             :     {
    2818    40469408 :       side->set_interior_parent(this);
    2819    38964080 :       side->inherit_data_from(*this);
    2820   284734697 :       for (auto n : side->node_index_range())
    2821   244265289 :         side->set_node(n, this->node_ptr(Subclass::side_nodes_map[i][n]));
    2822             :     }
    2823    41325596 : }
    2824             : 
    2825             : 
    2826             : 
    2827             : template <typename Subclass, typename Mapclass>
    2828             : inline
    2829             : void
    2830   247066302 : Elem::simple_side_ptr (std::unique_ptr<Elem> & side,
    2831             :                        const unsigned int i,
    2832             :                        ElemType sidetype)
    2833             : {
    2834     4831760 :   libmesh_assert_less (i, this->n_sides());
    2835             : 
    2836   247066302 :   if (!side.get() || side->type() != sidetype)
    2837             :     {
    2838       17462 :       Subclass & real_me = cast_ref<Subclass&>(*this);
    2839     1340692 :       side = real_me.Subclass::side_ptr(i);
    2840             :     }
    2841             :   else
    2842             :     {
    2843   251159115 :       side->subdomain_id() = this->subdomain_id();
    2844             : 
    2845   839007119 :       for (auto n : side->node_index_range())
    2846   592619894 :         side->set_node(n, this->node_ptr(Mapclass::side_nodes_map[i][n]));
    2847             :     }
    2848   247066302 : }
    2849             : 
    2850             : 
    2851             : 
    2852             : inline
    2853             : std::unique_ptr<const Elem>
    2854    48352376 : Elem::build_edge_ptr (const unsigned int i) const
    2855             : {
    2856             :   // Call the non-const version of this function, return the result as
    2857             :   // a std::unique_ptr<const Elem>.
    2858    28866390 :   Elem * me = const_cast<Elem *>(this);
    2859    50461417 :   return me->build_edge_ptr(i);
    2860             : }
    2861             : 
    2862             : 
    2863             : 
    2864             : inline
    2865             : void
    2866      151988 : Elem::build_edge_ptr (std::unique_ptr<const Elem> & elem,
    2867             :                       const unsigned int i) const
    2868             : {
    2869             :   // Hand off to the non-const version of this function
    2870        4508 :   Elem * me = const_cast<Elem *>(this);
    2871        9016 :   std::unique_ptr<Elem> e {const_cast<Elem *>(elem.release())};
    2872      151988 :   me->build_edge_ptr(e, i);
    2873      147480 :   elem = std::move(e);
    2874      151988 : }
    2875             : 
    2876             : 
    2877             : template <typename Edgeclass, typename Subclass>
    2878             : inline
    2879             : std::unique_ptr<Elem>
    2880    27678359 : Elem::simple_build_edge_ptr (const unsigned int i)
    2881             : {
    2882     9611726 :   libmesh_assert_less (i, this->n_edges());
    2883             : 
    2884    27678359 :   std::unique_ptr<Elem> edge = std::make_unique<Edgeclass>();
    2885             : 
    2886    97081094 :   for (auto n : edge->node_index_range())
    2887    69402735 :     edge->set_node(n, this->node_ptr(Subclass::edge_nodes_map[i][n]));
    2888             : 
    2889    27678359 :   edge->set_interior_parent(this);
    2890    25892105 :   edge->inherit_data_from(*this);
    2891             : 
    2892    27678359 :   return edge;
    2893           0 : }
    2894             : 
    2895             : 
    2896             : 
    2897             : 
    2898             : template <typename Subclass>
    2899             : inline
    2900             : void
    2901      182246 : Elem::simple_build_edge_ptr (std::unique_ptr<Elem> & edge,
    2902             :                              const unsigned int i,
    2903             :                              ElemType edgetype)
    2904             : {
    2905        6068 :   libmesh_assert_less (i, this->n_edges());
    2906             : 
    2907      182246 :   if (!edge.get() || edge->type() != edgetype)
    2908             :     {
    2909          35 :       Subclass & real_me = cast_ref<Subclass&>(*this);
    2910        1141 :       edge = real_me.Subclass::build_edge_ptr(i);
    2911             :     }
    2912             :   else
    2913             :     {
    2914      175625 :       edge->inherit_data_from(*this);
    2915      591520 :       for (auto n : edge->node_index_range())
    2916      409862 :         edge->set_node(n, this->node_ptr(Subclass::edge_nodes_map[i][n]));
    2917             :     }
    2918      182246 : }
    2919             : 
    2920             : 
    2921             : 
    2922             : inline
    2923          39 : bool Elem::on_boundary () const
    2924             : {
    2925             :   // By convention, the element is on the boundary
    2926             :   // if it has a nullptr neighbor.
    2927         351 :   return this->has_neighbor(nullptr);
    2928             : }
    2929             : 
    2930             : 
    2931             : 
    2932             : inline
    2933   155383942 : unsigned int Elem::which_neighbor_am_i (const Elem * e) const
    2934             : {
    2935     8647523 :   libmesh_assert(e);
    2936             : 
    2937     8647523 :   const Elem * eparent = e;
    2938             : 
    2939   156329120 :   while (eparent->level() > this->level())
    2940             :     {
    2941      339112 :       eparent = eparent->parent();
    2942      328992 :       libmesh_assert(eparent);
    2943             :     }
    2944             : 
    2945   408216876 :   for (auto s : make_range(this->n_sides()))
    2946   408209429 :     if (this->neighbor_ptr(s) == eparent)
    2947     8647523 :       return s;
    2948             : 
    2949           0 :   return libMesh::invalid_uint;
    2950             : }
    2951             : 
    2952             : 
    2953             : 
    2954             : inline
    2955   179824884 : bool Elem::active() const
    2956             : {
    2957             : #ifdef LIBMESH_ENABLE_AMR
    2958  4183419220 :   if ((this->refinement_flag() == INACTIVE) ||
    2959   123776941 :       (this->refinement_flag() == COARSEN_INACTIVE))
    2960    56174081 :     return false;
    2961             :   else
    2962   123650803 :     return true;
    2963             : #else
    2964             :   return true;
    2965             : #endif
    2966             : }
    2967             : 
    2968             : 
    2969             : 
    2970             : 
    2971             : 
    2972             : inline
    2973   963393160 : bool Elem::subactive() const
    2974             : {
    2975             : #ifdef LIBMESH_ENABLE_AMR
    2976    48595101 :   if (this->active())
    2977    35467202 :     return false;
    2978    13127899 :   if (!this->has_children())
    2979     3944800 :     return true;
    2980   285538348 :   for (const Elem * my_ancestor = this->parent();
    2981   820179410 :        my_ancestor != nullptr;
    2982    30752999 :        my_ancestor = my_ancestor->parent())
    2983    26663489 :     if (my_ancestor->active())
    2984       10516 :       return true;
    2985             : #endif
    2986             : 
    2987     9172583 :   return false;
    2988             : }
    2989             : 
    2990             : 
    2991             : 
    2992             : inline
    2993    49504673 : bool Elem::has_children() const
    2994             : {
    2995             : #ifdef LIBMESH_ENABLE_AMR
    2996   853661593 :   if (!_children)
    2997     9829108 :     return false;
    2998             :   else
    2999    39675565 :     return true;
    3000             : #else
    3001             :   return false;
    3002             : #endif
    3003             : }
    3004             : 
    3005             : 
    3006             : inline
    3007             : bool Elem::has_ancestor_children() const
    3008             : {
    3009             : #ifdef LIBMESH_ENABLE_AMR
    3010             :   if (!_children)
    3011             :     return false;
    3012             :   else
    3013             :     for (auto & c : child_ref_range())
    3014             :       if (c.has_children())
    3015             :         return true;
    3016             : #endif
    3017             :   return false;
    3018             : }
    3019             : 
    3020             : 
    3021             : 
    3022             : inline
    3023        1630 : bool Elem::is_ancestor_of(const Elem *
    3024             : #ifdef LIBMESH_ENABLE_AMR
    3025             :                           descendant
    3026             : #endif
    3027             :                           ) const
    3028             : {
    3029             : #ifdef LIBMESH_ENABLE_AMR
    3030        1630 :   const Elem * e = descendant;
    3031       12722 :   while (e)
    3032             :     {
    3033       10734 :       if (this == e)
    3034        1630 :         return true;
    3035         466 :       e = e->parent();
    3036             :     }
    3037             : #endif
    3038           0 :   return false;
    3039             : }
    3040             : 
    3041             : 
    3042             : 
    3043             : inline
    3044  1248289238 : const Elem * Elem::parent () const
    3045             : {
    3046  2626986650 :   return _elemlinks[0];
    3047             : }
    3048             : 
    3049             : 
    3050             : 
    3051             : inline
    3052    80082075 : Elem * Elem::parent ()
    3053             : {
    3054   974851154 :   return _elemlinks[0];
    3055             : }
    3056             : 
    3057             : 
    3058             : 
    3059             : inline
    3060      347654 : void Elem::set_parent (Elem * p)
    3061             : {
    3062             :   // We no longer support using parent() as interior_parent()
    3063      347654 :   libmesh_assert_equal_to(this->dim(), p ? p->dim() : this->dim());
    3064    10130201 :   _elemlinks[0] = p;
    3065      530509 : }
    3066             : 
    3067             : 
    3068             : 
    3069             : inline
    3070      801247 : const Elem * Elem::top_parent () const
    3071             : {
    3072      801247 :   const Elem * tp = this;
    3073             : 
    3074             :   // Keep getting the element's parent
    3075             :   // until that parent is at level-0
    3076     4304822 :   while (tp->parent() != nullptr)
    3077     1998084 :     tp = tp->parent();
    3078             : 
    3079      801247 :   libmesh_assert(tp);
    3080      801247 :   libmesh_assert_equal_to (tp->level(), 0);
    3081             : 
    3082      801247 :   return tp;
    3083             : }
    3084             : 
    3085             : 
    3086             : 
    3087             : inline
    3088  8048334485 : unsigned int Elem::level() const
    3089             : {
    3090             : #ifdef LIBMESH_ENABLE_AMR
    3091             : 
    3092             :   // if I don't have a parent I was
    3093             :   // created directly from file
    3094             :   // or by the user, so I am a
    3095             :   // level-0 element
    3096 22145099359 :   if (this->parent() == nullptr)
    3097   111258352 :     return 0;
    3098             : 
    3099             :   // if the parent and this element are of different
    3100             :   // dimensionality we are at the same level as
    3101             :   // the parent (e.g. we are the 2D side of a
    3102             :   // 3D element)
    3103 14419692893 :   if (this->dim() != this->parent()->dim())
    3104           0 :     return this->parent()->level();
    3105             : 
    3106             :   // otherwise we are at a level one
    3107             :   // higher than our parent
    3108 14552508216 :   return (this->parent()->level() + 1);
    3109             : 
    3110             : #else
    3111             : 
    3112             :   // Without AMR all elements are
    3113             :   // at level 0.
    3114             :   return 0;
    3115             : 
    3116             : #endif
    3117             : }
    3118             : 
    3119             : 
    3120             : 
    3121             : inline
    3122  2268217944 : unsigned int Elem::p_level() const
    3123             : {
    3124             : #ifdef LIBMESH_ENABLE_AMR
    3125 57822336316 :   return _p_level;
    3126             : #else
    3127             :   return 0;
    3128             : #endif
    3129             : }
    3130             : 
    3131             : 
    3132             : 
    3133             : inline
    3134   112433528 : ElemMappingType Elem::mapping_type () const
    3135             : {
    3136  1571554780 :   return static_cast<ElemMappingType>(_map_type);
    3137             : }
    3138             : 
    3139             : 
    3140             : 
    3141             : inline
    3142    48414308 : void Elem::set_mapping_type(const ElemMappingType type)
    3143             : {
    3144   311573510 :   _map_type = cast_int<unsigned char>(type);
    3145    48414308 : }
    3146             : 
    3147             : 
    3148             : 
    3149             : inline
    3150    34986160 : unsigned char Elem::mapping_data () const
    3151             : {
    3152   140811047 :   return _map_data;
    3153             : }
    3154             : 
    3155             : 
    3156             : 
    3157             : inline
    3158    48414301 : void Elem::set_mapping_data(const unsigned char data)
    3159             : {
    3160   311303195 :   _map_data = data;
    3161    48685977 : }
    3162             : 
    3163             : 
    3164             : 
    3165             : #ifdef LIBMESH_ENABLE_AMR
    3166             : 
    3167             : inline
    3168           0 : const Elem * Elem::raw_child_ptr (unsigned int i) const
    3169             : {
    3170        4240 :   if (!_children)
    3171           0 :     return nullptr;
    3172             : 
    3173        2880 :   return _children[i];
    3174             : }
    3175             : 
    3176             : inline
    3177    95094703 : const Elem * Elem::child_ptr (unsigned int i) const
    3178             : {
    3179    95094703 :   libmesh_assert(_children);
    3180    95094703 :   libmesh_assert(_children[i]);
    3181             : 
    3182   350458517 :   return _children[i];
    3183             : }
    3184             : 
    3185             : inline
    3186     2239934 : Elem * Elem::child_ptr (unsigned int i)
    3187             : {
    3188     2239934 :   libmesh_assert(_children);
    3189     2239934 :   libmesh_assert(_children[i]);
    3190             : 
    3191   329301889 :   return _children[i];
    3192             : }
    3193             : 
    3194             : 
    3195             : inline
    3196      134936 : void Elem::set_child (unsigned int c, Elem * elem)
    3197             : {
    3198      134936 :   libmesh_assert (this->has_children());
    3199             : 
    3200    52339099 :   _children[c] = elem;
    3201    27104089 : }
    3202             : 
    3203             : 
    3204             : 
    3205             : inline
    3206   112901817 : unsigned int Elem::which_child_am_i (const Elem * e) const
    3207             : {
    3208    29008140 :   libmesh_assert(e);
    3209    29008140 :   libmesh_assert (this->has_children());
    3210             : 
    3211   112901817 :   unsigned int nc = this->n_children();
    3212   300420790 :   for (unsigned int c=0; c != nc; c++)
    3213   300420790 :     if (this->child_ptr(c) == e)
    3214    29008140 :       return c;
    3215             : 
    3216           0 :   libmesh_error_msg("ERROR:  which_child_am_i() was called with a non-child!");
    3217             : 
    3218             :   return libMesh::invalid_uint;
    3219             : }
    3220             : 
    3221             : 
    3222             : 
    3223             : inline
    3224   355878374 : Elem::RefinementState Elem::refinement_flag () const
    3225             : {
    3226  5668517355 :   return static_cast<RefinementState>(_rflag);
    3227             : }
    3228             : 
    3229             : 
    3230             : 
    3231             : inline
    3232     3751870 : void Elem::set_refinement_flag(RefinementState rflag)
    3233             : {
    3234    62771328 :   _rflag = cast_int<unsigned char>(rflag);
    3235    13638499 : }
    3236             : 
    3237             : 
    3238             : 
    3239             : inline
    3240    96229096 : Elem::RefinementState Elem::p_refinement_flag () const
    3241             : {
    3242   270499932 :   return static_cast<RefinementState>(_pflag);
    3243             : }
    3244             : 
    3245             : 
    3246             : 
    3247             : inline
    3248     2591748 : void Elem::set_p_refinement_flag(RefinementState pflag)
    3249             : {
    3250     2591399 :   if (this->p_level() == 0)
    3251     2586135 :     libmesh_assert_not_equal_to
    3252             :       (pflag, Elem::JUST_REFINED);
    3253             : 
    3254    41344973 :   _pflag = cast_int<unsigned char>(pflag);
    3255    28190669 : }
    3256             : 
    3257             : 
    3258             : 
    3259             : inline
    3260           0 : unsigned int Elem::max_descendant_p_level () const
    3261             : {
    3262             :   // This is undefined for subactive elements,
    3263             :   // which have no active descendants
    3264           0 :   libmesh_assert (!this->subactive());
    3265           0 :   if (this->active())
    3266           0 :     return this->p_level();
    3267             : 
    3268           0 :   unsigned int max_p_level = _p_level;
    3269           0 :   for (auto & c : child_ref_range())
    3270           0 :     max_p_level = std::max(max_p_level,
    3271           0 :                            c.max_descendant_p_level());
    3272           0 :   return max_p_level;
    3273             : }
    3274             : 
    3275             : 
    3276             : 
    3277             : inline
    3278    48395890 : void Elem::hack_p_level(unsigned int p)
    3279             : {
    3280    48395890 :   if (p == 0)
    3281    48312479 :     libmesh_assert_not_equal_to
    3282             :       (this->p_refinement_flag(), Elem::JUST_REFINED);
    3283             : 
    3284   370111901 :   _p_level = cast_int<unsigned char>(p);
    3285   304848707 : }
    3286             : 
    3287             : 
    3288             : inline
    3289        4106 : void Elem::hack_p_level_and_refinement_flag (unsigned int p,
    3290             :                                              RefinementState pflag)
    3291             : {
    3292    38528085 :   _pflag = cast_int<unsigned char>(pflag);
    3293        4106 :   this->hack_p_level(p);
    3294    38508754 : }
    3295             : 
    3296             : #endif // ifdef LIBMESH_ENABLE_AMR
    3297             : 
    3298             : 
    3299             : inline
    3300      290819 : void Elem::orient(BoundaryInfo * boundary_info)
    3301             : {
    3302      304083 :   if (this->is_flipped())
    3303      138941 :     this->flip(boundary_info);
    3304      290819 : }
    3305             : 
    3306             : 
    3307             : inline
    3308       43994 : dof_id_type Elem::compute_key (dof_id_type n0)
    3309             : {
    3310       43994 :   return n0;
    3311             : }
    3312             : 
    3313             : 
    3314             : 
    3315             : inline
    3316     6788047 : dof_id_type Elem::compute_key (dof_id_type n0,
    3317             :                                dof_id_type n1)
    3318             : {
    3319             :   // Order the two so that n0 < n1
    3320   184618301 :   if (n0 > n1) std::swap (n0, n1);
    3321             : 
    3322   184618301 :   return Utility::hashword2(n0, n1);
    3323             : }
    3324             : 
    3325             : 
    3326             : 
    3327             : inline
    3328    55015624 : dof_id_type Elem::compute_key (dof_id_type n0,
    3329             :                                dof_id_type n1,
    3330             :                                dof_id_type n2)
    3331             : {
    3332    55015624 :   std::array<dof_id_type, 3> array = {{n0, n1, n2}};
    3333     1795464 :   std::sort(array.begin(), array.end());
    3334    56811088 :   return Utility::hashword(array);
    3335             : }
    3336             : 
    3337             : 
    3338             : 
    3339             : inline
    3340    36497379 : dof_id_type Elem::compute_key (dof_id_type n0,
    3341             :                                dof_id_type n1,
    3342             :                                dof_id_type n2,
    3343             :                                dof_id_type n3)
    3344             : {
    3345    36497379 :   std::array<dof_id_type, 4> array = {{n0, n1, n2, n3}};
    3346     1025658 :   std::sort(array.begin(), array.end());
    3347    37523037 :   return Utility::hashword(array);
    3348             : }
    3349             : 
    3350             : 
    3351             : 
    3352             : inline
    3353   235707797 : void Elem::inherit_data_from (const Elem & src)
    3354             : {
    3355    62599688 :   this->set_mapping_type(src.mapping_type());
    3356    62599688 :   this->set_mapping_data(src.mapping_data());
    3357   251182045 :   this->subdomain_id() = src.subdomain_id();
    3358    62599688 :   this->processor_id(src.processor_id());
    3359             : #ifdef LIBMESH_ENABLE_AMR
    3360   251182045 :   this->set_p_level(src.p_level());
    3361             : #endif
    3362   235707797 : }
    3363             : 
    3364             : 
    3365             : 
    3366             : /**
    3367             :  * The definition of the protected nested SideIter class.
    3368             :  */
    3369           0 : class Elem::SideIter
    3370             : {
    3371             : public:
    3372             :   // Constructor with arguments.
    3373           0 :   SideIter(const unsigned int side_number,
    3374             :            Elem * parent)
    3375           0 :     : _side(),
    3376           0 :       _side_ptr(nullptr),
    3377           0 :       _parent(parent),
    3378           0 :       _side_number(side_number)
    3379           0 :   {}
    3380             : 
    3381             : 
    3382             :   // Empty constructor.
    3383             :   SideIter()
    3384             :     : _side(),
    3385             :       _side_ptr(nullptr),
    3386             :       _parent(nullptr),
    3387             :       _side_number(libMesh::invalid_uint)
    3388             :   {}
    3389             : 
    3390             : 
    3391             :   // Copy constructor
    3392           0 :   SideIter(const SideIter & other)
    3393           0 :     : _side(),
    3394           0 :       _side_ptr(nullptr),
    3395           0 :       _parent(other._parent),
    3396           0 :       _side_number(other._side_number)
    3397           0 :   {}
    3398             : 
    3399             : 
    3400             :   // op=
    3401             :   SideIter & operator=(const SideIter & other)
    3402             :   {
    3403             :     this->_parent      = other._parent;
    3404             :     this->_side_number = other._side_number;
    3405             :     return *this;
    3406             :   }
    3407             : 
    3408             :   // unary op*
    3409           0 :   Elem *& operator*() const
    3410             :   {
    3411             :     // Set the std::unique_ptr
    3412           0 :     this->_update_side_ptr();
    3413             : 
    3414             :     // Return a reference to _side_ptr
    3415           0 :     return this->_side_ptr;
    3416             :   }
    3417             : 
    3418             :   // op++
    3419           0 :   SideIter & operator++()
    3420             :   {
    3421           0 :     ++_side_number;
    3422           0 :     return *this;
    3423             :   }
    3424             : 
    3425             :   // op==  Two side iterators are equal if they have
    3426             :   // the same side number and the same parent element.
    3427           0 :   bool operator == (const SideIter & other) const
    3428             :   {
    3429           0 :     return (this->_side_number == other._side_number &&
    3430           0 :             this->_parent      == other._parent);
    3431             :   }
    3432             : 
    3433             : 
    3434             :   // Consults the parent Elem to determine if the side
    3435             :   // is a boundary side.  Note: currently side N is a
    3436             :   // boundary side if neighbor N is nullptr.  Be careful,
    3437             :   // this could possibly change in the future?
    3438           0 :   bool side_on_boundary() const
    3439             :   {
    3440           0 :     return this->_parent->neighbor_ptr(_side_number) == nullptr;
    3441             :   }
    3442             : 
    3443             : private:
    3444             :   // Update the _side pointer by building the correct side.
    3445             :   // This has to be called before dereferencing.
    3446           0 :   void _update_side_ptr() const
    3447             :   {
    3448             :     // Construct new side, store in std::unique_ptr
    3449           0 :     this->_side = this->_parent->build_side_ptr(this->_side_number);
    3450             : 
    3451             :     // Also set our internal naked pointer.  Memory is still owned
    3452             :     // by the std::unique_ptr.
    3453           0 :     this->_side_ptr = _side.get();
    3454           0 :   }
    3455             : 
    3456             :   // std::unique_ptr to the actual side, handles memory management for
    3457             :   // the sides which are created during the course of iteration.
    3458             :   mutable std::unique_ptr<Elem> _side;
    3459             : 
    3460             :   // Raw pointer needed to facilitate passing back to the user a
    3461             :   // reference to a non-temporary raw pointer in order to conform to
    3462             :   // the variant_filter_iterator interface.  It points to the same
    3463             :   // thing the std::unique_ptr "_side" above holds.  What happens if the user
    3464             :   // calls delete on the pointer passed back?  Well, this is an issue
    3465             :   // which is not addressed by the iterators in libMesh.  Basically it
    3466             :   // is a bad idea to ever call delete on an iterator from the library.
    3467             :   mutable Elem * _side_ptr;
    3468             : 
    3469             :   // Pointer to the parent Elem class which generated this iterator
    3470             :   Elem * _parent;
    3471             : 
    3472             :   // A counter variable which keeps track of the side number
    3473             :   unsigned int _side_number;
    3474             : };
    3475             : 
    3476             : 
    3477             : 
    3478             : 
    3479             : 
    3480             : 
    3481             : // Private implementation functions in the Elem class for the side iterators.
    3482             : // They have to come after the definition of the SideIter class.
    3483             : inline
    3484           0 : Elem::SideIter Elem::_first_side()
    3485             : {
    3486           0 :   return SideIter(0, this);
    3487             : }
    3488             : 
    3489             : 
    3490             : 
    3491             : inline
    3492           0 : Elem::SideIter Elem::_last_side()
    3493             : {
    3494           0 :   return SideIter(this->n_neighbors(), this);
    3495             : }
    3496             : 
    3497             : 
    3498             : 
    3499             : 
    3500             : /**
    3501             :  * The definition of the struct used for iterating over sides.
    3502             :  */
    3503             : struct
    3504             : Elem::side_iterator : variant_filter_iterator<Elem::Predicate, Elem *>
    3505             : {
    3506             :   // Templated forwarding ctor -- forwards to appropriate variant_filter_iterator ctor
    3507             :   template <typename PredType, typename IterType>
    3508           0 :   side_iterator (const IterType & d,
    3509             :                  const IterType & e,
    3510             :                  const PredType & p ) :
    3511           0 :     variant_filter_iterator<Elem::Predicate, Elem *>(d,e,p) {}
    3512             : };
    3513             : 
    3514             : 
    3515             : 
    3516             : inline
    3517    86944614 : SimpleRange<Elem::NeighborPtrIter> Elem::neighbor_ptr_range()
    3518             : {
    3519    90236036 :   return {_elemlinks+1, _elemlinks + 1 + this->n_neighbors()};
    3520             : }
    3521             : 
    3522             : 
    3523             : inline
    3524   403258739 : SimpleRange<Elem::ConstNeighborPtrIter> Elem::neighbor_ptr_range() const
    3525             : {
    3526   404811762 :   return {_elemlinks+1, _elemlinks + 1 + this->n_neighbors()};
    3527             : }
    3528             : 
    3529             : } // namespace libMesh
    3530             : 
    3531             : 
    3532             : // Helper function for default caches in Elem subclasses
    3533             : 
    3534             : #define LIBMESH_ENABLE_TOPOLOGY_CACHES                                  \
    3535             :   virtual                                                               \
    3536             :   std::vector<std::vector<std::vector<std::vector<std::pair<unsigned char, unsigned char>>>>> & \
    3537             :   _get_bracketing_node_cache() const override                   \
    3538             :   {                                                                     \
    3539             :     static std::vector<std::vector<std::vector<std::vector<std::pair<unsigned char, unsigned char>>>>> c; \
    3540             :     return c;                                                           \
    3541             :   }                                                                     \
    3542             :                                                                         \
    3543             :   virtual                                                               \
    3544             :   std::vector<std::vector<std::vector<signed char>>> &                  \
    3545             :   _get_parent_indices_cache() const override                    \
    3546             :   {                                                                     \
    3547             :     static std::vector<std::vector<std::vector<signed char>>> c;        \
    3548             :     return c;                                                           \
    3549             :   }
    3550             : 
    3551             : 
    3552             : 
    3553             : 
    3554             : 
    3555             : 
    3556             : #endif // LIBMESH_ELEM_H

Generated by: LCOV version 1.14