LCOV - code coverage report
Current view: top level - include/geom - elem.h (source / functions) Hit Total Coverage
Test: libMesh/libmesh: #4506 (f0a4b5) with base 5112d2 Lines: 355 436 81.4 %
Date: 2026-07-29 17:31:56 Functions: 142 182 78.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14