LCOV - code coverage report
Current view: top level - include/geom - elem.h (source / functions) Hit Total Coverage
Test: libMesh/libmesh: #4250 (d9a2b9) with base 332933 Lines: 353 437 80.8 %
Date: 2025-09-29 21:21:55 Functions: 140 182 76.9 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14