elem.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#ifndef LIBMESH_ELEM_H
#define LIBMESH_ELEM_H

// Local includes
#include "libmesh/libmesh_common.h"
#include "libmesh/bounding_box.h"
#include "libmesh/dof_object.h"
#include "libmesh/id_types.h"
#include "libmesh/reference_counted_object.h"
#include "libmesh/node.h"
#include "libmesh/enum_elem_type.h" // INVALID_ELEM
#include "libmesh/multi_predicates.h"
#include "libmesh/pointer_to_pointer_iter.h"
#include "libmesh/int_range.h"
#include "libmesh/simple_range.h"
#include "libmesh/variant_filter_iterator.h"
#include "libmesh/hashword.h" // Used in compute_key() functions

// C++ includes
#include <algorithm>
#include <cstddef>
#include <iostream>
#include <limits.h> // CHAR_BIT
#include <set>
#include <vector>
#include <memory>
#include <array>

namespace libMesh
{

// Forward declarations
class BoundaryInfo;
class Elem;
class MeshBase;
class MeshRefinement;
#ifdef LIBMESH_ENABLE_PERIODIC
class PeriodicBoundaries;
class PointLocatorBase;
#endif
template <class SideType, class ParentType>
class Side;
enum ElemQuality : int;
enum IOPackage : int;
enum Order : int;


class Elem : public ReferenceCountedObject<Elem>,
             public DofObject
{
protected:

  Elem (const unsigned int n_nodes,
        const unsigned int n_sides,
        Elem * parent,
        Elem ** elemlinkdata,
        Node ** nodelinkdata);

public:

  Elem (Elem &&) = delete;
  Elem (const Elem &) = delete;
  Elem & operator= (const Elem &) = delete;
  Elem & operator= (Elem &&) = delete;

  virtual ~Elem() = default;

  const Point & point (const unsigned int i) const;

  Point & point (const unsigned int i);

  virtual Point master_point (const unsigned int i) const = 0;

  dof_id_type node_id (const unsigned int i) const;

  unsigned int local_node (const dof_id_type i) const;

  unsigned int get_node_index (const Node * node_ptr) const;

  const Node * const * get_nodes () const;

  const Node * node_ptr (const unsigned int i) const;

  Node * node_ptr (const unsigned int i);

  const Node & node_ref (const unsigned int i) const;

  Node & node_ref (const unsigned int i);

#ifdef LIBMESH_ENABLE_DEPRECATED

  virtual Node * & set_node (const unsigned int i);
#endif // LIBMESH_ENABLE_DEPRECATED

  virtual void set_node (const unsigned int i,
                         Node * node);

  class NodeRefIter;
  class ConstNodeRefIter;

  SimpleRange<NodeRefIter> node_ref_range();

  SimpleRange<ConstNodeRefIter> node_ref_range() const;

  subdomain_id_type subdomain_id () const;

  subdomain_id_type & subdomain_id ();

  static constexpr subdomain_id_type invalid_subdomain_id
    = std::numeric_limits<subdomain_id_type>::max();

  virtual bool runtime_topology() const { return false; }

  const Elem * reference_elem () const;

  virtual dof_id_type key (const unsigned int s) const = 0;

  virtual dof_id_type low_order_key (const unsigned int s) const = 0;

  virtual dof_id_type key () const;

  bool operator == (const Elem & rhs) const;

  bool operator != (const Elem & rhs) const;

  bool topologically_equal (const Elem & rhs) const;

  const Elem * neighbor_ptr (unsigned int i) const;

  Elem * neighbor_ptr (unsigned int i);

  typedef Elem * const *       NeighborPtrIter;
  typedef const Elem * const * ConstNeighborPtrIter;

  SimpleRange<NeighborPtrIter> neighbor_ptr_range();

  SimpleRange<ConstNeighborPtrIter> neighbor_ptr_range() const;

#ifdef LIBMESH_ENABLE_PERIODIC

  const Elem * topological_neighbor (const unsigned int i,
                                     const MeshBase & mesh,
                                     const PointLocatorBase & point_locator,
                                     const PeriodicBoundaries * pb) const;

  Elem * topological_neighbor (const unsigned int i,
                               MeshBase & mesh,
                               const PointLocatorBase & point_locator,
                               const PeriodicBoundaries * pb);

  bool has_topological_neighbor (const Elem * elem,
                                 const MeshBase & mesh,
                                 const PointLocatorBase & point_locator,
                                 const PeriodicBoundaries * pb) const;
#endif

  void set_neighbor (const unsigned int i, Elem * n);

  bool has_neighbor (const Elem * elem) const;

  Elem * child_neighbor (Elem * elem);

  const Elem * child_neighbor (const Elem * elem) const;

  bool on_boundary () const;

  bool is_semilocal (const processor_id_type my_pid) const;

  unsigned int which_neighbor_am_i(const Elem * e) const;

  unsigned int which_side_am_i(const Elem * e) const;

  virtual unsigned int local_side_node(unsigned int side,
                                       unsigned int side_node) const = 0;

  virtual unsigned int local_edge_node(unsigned int edge,
                                       unsigned int edge_node) const = 0;

  bool contains_vertex_of(const Elem * e, bool mesh_connection=false) const;

  bool contains_edge_of(const Elem * e) const;

  void find_point_neighbors(const Point & p,
                            std::set<const Elem *> & neighbor_set) const;

  void find_point_neighbors(std::set<const Elem *> & neighbor_set) const;

  void find_point_neighbors(std::set<const Elem *> & neighbor_set,
                            const Elem * start_elem) const;

  void find_point_neighbors(std::set<Elem *> & neighbor_set,
                            Elem * start_elem);

  void find_edge_neighbors(const Point & p1,
                           const Point & p2,
                           std::set<const Elem *> & neighbor_set) const;

  void find_edge_neighbors(std::set<const Elem *> & neighbor_set) const;

  void find_interior_neighbors(std::set<const Elem *> & neighbor_set) const;

  void find_interior_neighbors(std::set<Elem *> & neighbor_set);

  void remove_links_to_me ();

  void make_links_to_me_remote ();

  void make_links_to_me_local (unsigned int n, unsigned int neighbor_side);

  virtual bool is_remote () const
  { return false; }

  virtual void connectivity(const unsigned int sc,
                            const IOPackage iop,
                            std::vector<dof_id_type> & conn) const = 0;

  void write_connectivity (std::ostream & out,
                           const IOPackage iop) const;

  virtual ElemType type () const = 0;

  static const unsigned int type_to_dim_map[INVALID_ELEM];

  virtual unsigned short dim () const = 0;

  static const unsigned int type_to_n_nodes_map[INVALID_ELEM];

  virtual unsigned int n_nodes () const = 0;

  static const unsigned int max_n_nodes = 27;

  IntRange<unsigned short> node_index_range () const;

  virtual unsigned int n_nodes_in_child (unsigned int /*c*/) const
  { return this->n_nodes(); }

  static const unsigned int type_to_n_sides_map[INVALID_ELEM];

  virtual unsigned int n_sides () const = 0;

  virtual ElemType side_type (const unsigned int s) const = 0;

  virtual Point side_vertex_average_normal(const unsigned int s) const;

  IntRange<unsigned short> side_index_range () const;

  unsigned int n_neighbors () const
  { return this->n_sides(); }

  virtual unsigned int n_vertices () const = 0;

  virtual unsigned int n_edges () const = 0;

  IntRange<unsigned short> edge_index_range () const;

  static const unsigned int type_to_n_edges_map[INVALID_ELEM];

  virtual unsigned int n_faces () const = 0;

  IntRange<unsigned short> face_index_range () const;

  virtual unsigned int n_children () const = 0;

  virtual bool is_vertex(const unsigned int i) const = 0;

  virtual bool is_vertex_on_child (unsigned int /*c*/,
                                   unsigned int n) const
  { return this->is_vertex(n); }

  virtual bool is_vertex_on_parent(unsigned int c,
                                   unsigned int n) const;

  virtual bool is_edge(const unsigned int i) const = 0;

  virtual bool is_face(const unsigned int i) const = 0;

  bool is_internal(const unsigned int i) const;

  virtual bool is_node_on_side(const unsigned int n,
                               const unsigned int s) const = 0;

  virtual std::vector<unsigned int> nodes_on_side(const unsigned int /*s*/) const = 0;

  virtual std::vector<unsigned int> nodes_on_edge(const unsigned int /*e*/) const = 0;

  virtual std::vector<unsigned int> sides_on_edge(const unsigned int /*e*/) const = 0;

  virtual std::vector<unsigned int> edges_adjacent_to_node(const unsigned int /*n*/) const = 0;

  virtual bool is_node_on_edge(const unsigned int n,
                               const unsigned int e) const = 0;

  virtual bool is_edge_on_side(const unsigned int e,
                               const unsigned int s) const = 0;

  virtual unsigned int opposite_side(const unsigned int s) const;

  virtual unsigned int opposite_node(const unsigned int n,
                                     const unsigned int s) const;

  virtual unsigned int n_sub_elem () const = 0;

  virtual std::unique_ptr<Elem> side_ptr (unsigned int i) = 0;
  std::unique_ptr<const Elem> side_ptr (unsigned int i) const;

  virtual void side_ptr (std::unique_ptr<Elem> & side, const unsigned int i) = 0;
  void side_ptr (std::unique_ptr<const Elem> & side, const unsigned int i) const;

  virtual std::unique_ptr<Elem> build_side_ptr (const unsigned int i) = 0;
  std::unique_ptr<const Elem> build_side_ptr (const unsigned int i) const;

#ifdef LIBMESH_ENABLE_DEPRECATED
  /*
   * Older versions of libMesh supported a "proxy" option here.
   */
  virtual std::unique_ptr<Elem> build_side_ptr (const unsigned int i, bool proxy)
  { if (proxy) libmesh_error(); libmesh_deprecated(); return this->build_side_ptr(i); }

  std::unique_ptr<const Elem> build_side_ptr (const unsigned int i, bool proxy) const
  { if (proxy) libmesh_error(); libmesh_deprecated(); return this->build_side_ptr(i); }
#endif

  virtual void build_side_ptr (std::unique_ptr<Elem> & side, const unsigned int i) = 0;
  void build_side_ptr (std::unique_ptr<const Elem> & side, const unsigned int i) const;

  virtual std::unique_ptr<Elem> build_edge_ptr (const unsigned int i) = 0;
  std::unique_ptr<const Elem> build_edge_ptr (const unsigned int i) const;

  virtual void build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i) = 0;
  void build_edge_ptr (std::unique_ptr<const Elem> & edge, const unsigned int i) const;

  static const Order type_to_default_order_map[INVALID_ELEM];

  virtual Order default_order () const = 0;

  virtual Order supported_nodal_order() const { return default_order(); }

  virtual Order default_side_order () const { return default_order(); }

  virtual Point true_centroid () const;

  Point vertex_average () const;

  virtual Point quasicircumcenter () const
  { libmesh_not_implemented(); }

  virtual Real hmin () const;

  virtual Real hmax () const;

  virtual Real volume () const;

  virtual BoundingBox loose_bounding_box () const;

  virtual Real quality (const ElemQuality q) const;

  virtual std::pair<Real,Real> qual_bounds (const ElemQuality) const
  { libmesh_not_implemented(); return std::make_pair(0.,0.); }

  virtual bool contains_point (const Point & p, Real tol=TOLERANCE) const;

  virtual bool on_reference_element(const Point & p,
                                    const Real eps = TOLERANCE) const = 0;

  virtual bool close_to_point(const Point & p, Real tol) const;

  bool positive_edge_orientation(const unsigned int i) const;

  bool positive_face_orientation(const unsigned int i) const;

  bool relative_edge_face_order(const unsigned int e, const unsigned int s) const;

  void inherit_data_from(const Elem & src);

private:
  bool point_test(const Point & p, Real box_tol, Real map_tol) const;

public:
  virtual bool has_affine_map () const { return false; }

  virtual bool has_invertible_map(Real tol = TOLERANCE*TOLERANCE) const;

  virtual bool is_linear () const { return false; }

  void print_info (std::ostream & os=libMesh::out) const;

  std::string get_info () const;

  bool active () const;

  bool ancestor () const;

  bool subactive () const;

  bool has_children () const;

  bool has_ancestor_children () const;

  bool is_ancestor_of(const Elem * descendant) const;

  const Elem * parent () const;

  Elem * parent ();

  void set_parent (Elem * p);

  const Elem * top_parent () const;

  const Elem * interior_parent () const;

  Elem * interior_parent ();

  void set_interior_parent (Elem * p);

  Real length (const unsigned int n1,
               const unsigned int n2) const;

  virtual unsigned int n_second_order_adjacent_vertices (const unsigned int n) const;

  virtual unsigned short int second_order_adjacent_vertex (const unsigned int n,
                                                           const unsigned int v) const;

  virtual std::pair<unsigned short int, unsigned short int>
  second_order_child_vertex (const unsigned int n) const;

  static ElemType second_order_equivalent_type (const ElemType et,
                                                const bool full_ordered=true);

  static ElemType first_order_equivalent_type (const ElemType et);

  static ElemType complete_order_equivalent_type (const ElemType et);

  unsigned int level () const;

  unsigned int p_level () const;

  virtual bool is_child_on_side(const unsigned int c,
                                const unsigned int s) const = 0;

  ElemMappingType mapping_type () const;

  void set_mapping_type (const ElemMappingType type);

  unsigned char mapping_data () const;

  void set_mapping_data (const unsigned char data);


#ifdef LIBMESH_ENABLE_AMR

  enum RefinementState { COARSEN = 0,
                         DO_NOTHING,
                         REFINE,
                         JUST_REFINED,
                         JUST_COARSENED,
                         INACTIVE,
                         COARSEN_INACTIVE,
                         INVALID_REFINEMENTSTATE };

  const Elem * raw_child_ptr (unsigned int i) const;

  const Elem * child_ptr (unsigned int i) const;

  Elem * child_ptr (unsigned int i);

  class ChildRefIter;
  class ConstChildRefIter;

  SimpleRange<ChildRefIter> child_ref_range();

  SimpleRange<ConstChildRefIter> child_ref_range() const;

private:
  void set_child (unsigned int c, Elem * elem);

public:
  unsigned int which_child_am_i(const Elem * e) const;

  virtual bool is_child_on_edge(const unsigned int c,
                                const unsigned int e) const;

  void add_child (Elem * elem);

  void add_child (Elem * elem, unsigned int c);

  void replace_child (Elem * elem, unsigned int c);

  void family_tree (std::vector<const Elem *> & family,
                    bool reset = true) const;

  void family_tree (std::vector<Elem *> & family,
                    bool reset = true);

  void total_family_tree (std::vector<const Elem *> & family,
                          bool reset = true) const;

  void total_family_tree (std::vector<Elem *> & family,
                          bool reset = true);

  void active_family_tree (std::vector<const Elem *> & active_family,
                           bool reset = true) const;

  void active_family_tree (std::vector<Elem *> & active_family,
                           bool reset = true);

  void family_tree_by_side (std::vector<const Elem *> & family,
                            unsigned int side,
                            bool reset = true) const;

  void family_tree_by_side (std::vector<Elem *> & family,
                            unsigned int side,
                            bool reset = true);

  void total_family_tree_by_side (std::vector<const Elem *> & family,
                                  unsigned int side,
                                  bool reset = true) const;

  void total_family_tree_by_side (std::vector<Elem *> & family,
                                  unsigned int side,
                                  bool reset = true);

  void active_family_tree_by_side (std::vector<const Elem *> & family,
                                   unsigned int side,
                                   bool reset = true) const;

  void active_family_tree_by_side (std::vector<Elem *> & family,
                                   unsigned int side,
                                   bool reset = true);

  void family_tree_by_neighbor (std::vector<const Elem *> & family,
                                const Elem * neighbor,
                                bool reset = true) const;

  void family_tree_by_neighbor (std::vector<Elem *> & family,
                                Elem * neighbor,
                                bool reset = true);

  void total_family_tree_by_neighbor (std::vector<const Elem *> & family,
                                      const Elem * neighbor,
                                      bool reset = true) const;

  void total_family_tree_by_neighbor (std::vector<Elem *> & family,
                                      Elem * neighbor,
                                      bool reset = true);

  void family_tree_by_subneighbor (std::vector<const Elem *> & family,
                                   const Elem * neighbor,
                                   const Elem * subneighbor,
                                   bool reset = true) const;

  void family_tree_by_subneighbor (std::vector<Elem *> & family,
                                   Elem * neighbor,
                                   Elem * subneighbor,
                                   bool reset = true);

  void total_family_tree_by_subneighbor (std::vector<const Elem *> & family,
                                         const Elem * neighbor,
                                         const Elem * subneighbor,
                                         bool reset = true) const;

  void total_family_tree_by_subneighbor (std::vector<Elem *> & family,
                                         Elem * neighbor,
                                         Elem * subneighbor,
                                         bool reset = true);

  void active_family_tree_by_neighbor (std::vector<const Elem *> & family,
                                       const Elem * neighbor,
                                       bool reset = true) const;

  void active_family_tree_by_neighbor (std::vector<Elem *> & family,
                                       Elem * neighbor,
                                       bool reset = true);

  void active_family_tree_by_topological_neighbor (std::vector<const Elem *> & family,
                                                   const Elem * neighbor,
                                                   const MeshBase & mesh,
                                                   const PointLocatorBase & point_locator,
                                                   const PeriodicBoundaries * pb,
                                                   bool reset = true) const;

  void active_family_tree_by_topological_neighbor (std::vector<Elem *> & family,
                                                   Elem * neighbor,
                                                   const MeshBase & mesh,
                                                   const PointLocatorBase & point_locator,
                                                   const PeriodicBoundaries * pb,
                                                   bool reset = true);

  RefinementState refinement_flag () const;

  void set_refinement_flag (const RefinementState rflag);

  RefinementState p_refinement_flag () const;

  void set_p_refinement_flag (const RefinementState pflag);

  unsigned int max_descendant_p_level () const;

  unsigned int min_p_level_by_neighbor (const Elem * neighbor,
                                        unsigned int current_min) const;

  unsigned int min_new_p_level_by_neighbor (const Elem * neighbor,
                                            unsigned int current_min) const;

  void set_p_level (const unsigned int p);

  void hack_p_level (const unsigned int p);

  void hack_p_level_and_refinement_flag (const unsigned int p,
                                         RefinementState pflag);

  virtual void refine (MeshRefinement & mesh_refinement);

  void coarsen ();

  void contract ();

#endif

#ifndef NDEBUG

  void libmesh_assert_valid_neighbors() const;

  void libmesh_assert_valid_node_pointers() const;
#endif // !NDEBUG

  virtual unsigned int local_singular_node(const Point & /* p */, const Real /* tol */ = TOLERANCE*TOLERANCE) const
  { return invalid_uint; }

  virtual bool is_singular_node(unsigned int /* node_i */) const { return false; }

  virtual unsigned int center_node_on_side(const unsigned short side) const;

protected:

  class SideIter;

public:
  typedef Predicates::multi_predicate Predicate;

  struct side_iterator;

  side_iterator boundary_sides_begin();
  side_iterator boundary_sides_end();

private:
  SideIter _first_side();
  SideIter _last_side();

public:

#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS

  virtual bool infinite () const = 0;

  virtual bool is_mid_infinite_edge_node(const unsigned int /* n */) const
  { libmesh_assert (!this->infinite()); return false; }

  virtual Point origin () const { libmesh_not_implemented(); return Point(); }

#else

  static constexpr bool infinite () { return false; }

#endif

  static std::unique_ptr<Elem> build (const ElemType type,
                                      Elem * p=nullptr);

  static std::unique_ptr<Elem> build_with_id (const ElemType type,
                                              dof_id_type id);

  virtual std::unique_ptr<Elem> disconnected_clone () const;

  virtual unsigned int n_permutations() const = 0;

  virtual void permute(unsigned int perm_num) = 0;

  virtual void flip(BoundaryInfo * boundary_info) = 0;

   virtual bool is_flipped() const = 0;

  void orient(BoundaryInfo * boundary_info);

#ifdef LIBMESH_ENABLE_AMR

  virtual unsigned int as_parent_node (unsigned int c,
                                       unsigned int n) const;

  virtual
  const std::vector<std::pair<unsigned char, unsigned char>> &
  parent_bracketing_nodes(unsigned int c,
                          unsigned int n) const;

  virtual
  const std::vector<std::pair<dof_id_type, dof_id_type>>
  bracketing_nodes(unsigned int c,
                   unsigned int n) const;


  virtual Real embedding_matrix (const unsigned int child_num,
                                 const unsigned int child_node_num,
                                 const unsigned int parent_node_num) const = 0;

  virtual unsigned int embedding_matrix_version () const { return 0; }

#endif // LIBMESH_ENABLE_AMR


protected:

  static constexpr Real affine_tol = TOLERANCE*TOLERANCE;

  static dof_id_type compute_key (dof_id_type n0);

  static dof_id_type compute_key (dof_id_type n0,
                                  dof_id_type n1);

  static dof_id_type compute_key (dof_id_type n0,
                                  dof_id_type n1,
                                  dof_id_type n2);

  static dof_id_type compute_key (dof_id_type n0,
                                  dof_id_type n1,
                                  dof_id_type n2,
                                  dof_id_type n3);

  void swap2nodes(unsigned int n1, unsigned int n2)
  {
    Node * temp = this->node_ptr(n1);
    this->set_node(n1, this->node_ptr(n2));
    this->set_node(n2, temp);
  }

  void swap2neighbors(unsigned int n1, unsigned int n2)
  {
    Elem * temp = this->neighbor_ptr(n1);
    this->set_neighbor(n1, this->neighbor_ptr(n2));
    this->set_neighbor(n2, temp);
  }

  void swap2boundarysides(unsigned short s1, unsigned short s2,
                          BoundaryInfo * boundary_info) const;

  void swap2boundaryedges(unsigned short e1, unsigned short e2,
                          BoundaryInfo * boundary_info) const;

  void swap3nodes(unsigned int n1, unsigned int n2, unsigned int n3)
  {
    swap2nodes(n1, n2);
    swap2nodes(n2, n3);
  }

  void swap3neighbors(unsigned int n1, unsigned int n2,
                      unsigned int n3)
  {
    swap2neighbors(n1, n2);
    swap2neighbors(n2, n3);
  }

  void swap4nodes(unsigned int n1, unsigned int n2, unsigned int n3,
                  unsigned int n4)
  {
    swap3nodes(n1, n2, n3);
    swap2nodes(n3, n4);
  }

  void swap4neighbors(unsigned int n1, unsigned int n2,
                      unsigned int n3, unsigned int n4)
  {
    swap3neighbors(n1, n2, n3);
    swap2neighbors(n3, n4);
  }


  template <typename Sideclass, typename Subclass>
  std::unique_ptr<Elem>
  simple_build_side_ptr(const unsigned int i);

  template <typename Subclass>
  void simple_build_side_ptr(std::unique_ptr<Elem> & side,
                             const unsigned int i,
                             ElemType sidetype);

  template <typename Subclass, typename Mapclass>
  void simple_side_ptr(std::unique_ptr<Elem> & side,
                       const unsigned int i,
                       ElemType sidetype);

  template <typename Edgeclass, typename Subclass>
  std::unique_ptr<Elem>
  simple_build_edge_ptr(const unsigned int i);

  template <typename Subclass>
  void simple_build_edge_ptr(std::unique_ptr<Elem> & edge,
                             const unsigned int i,
                             ElemType edgetype);


#ifdef LIBMESH_ENABLE_AMR

  virtual
  std::vector<std::vector<std::vector<std::vector<std::pair<unsigned char, unsigned char>>>>> &
  _get_bracketing_node_cache() const
  {
    static std::vector<std::vector<std::vector<std::vector<std::pair<unsigned char, unsigned char>>>>> c;
    libmesh_error();
    return c;
  }

  virtual
  std::vector<std::vector<std::vector<signed char>>> &
  _get_parent_indices_cache() const
  {
    static std::vector<std::vector<std::vector<signed char>>> c;
    libmesh_error();
    return c;
  }

#endif // LIBMESH_ENABLE_AMR

public:

  void nullify_neighbors ();

protected:

  Node ** _nodes;

  Elem ** _elemlinks;

#ifdef LIBMESH_ENABLE_AMR

  std::unique_ptr<Elem *[]> _children;
#endif

  subdomain_id_type _sbd_id;

#ifdef LIBMESH_ENABLE_AMR

  unsigned char _rflag;

  unsigned char _pflag;

  unsigned char _p_level;
#endif

  unsigned char _map_type;

  unsigned char _map_data;
};



// ------------------------------------------------------------
// Elem helper classes
//
class
Elem::NodeRefIter : public PointerToPointerIter<Node>
{
public:
  NodeRefIter (Node * const * nodepp) : PointerToPointerIter<Node>(nodepp) {}
};


class
Elem::ConstNodeRefIter : public PointerToPointerIter<const Node>
{
public:
  ConstNodeRefIter (const Node * const * nodepp) : PointerToPointerIter<const Node>(nodepp) {}
};


#ifdef LIBMESH_ENABLE_AMR
class
Elem::ChildRefIter : public PointerToPointerIter<Elem>
{
public:
  ChildRefIter (Elem * const * childpp) : PointerToPointerIter<Elem>(childpp) {}
};


class
Elem::ConstChildRefIter : public PointerToPointerIter<const Elem>
{
public:
  ConstChildRefIter (const Elem * const * childpp) : PointerToPointerIter<const Elem>(childpp) {}
};



inline
SimpleRange<Elem::ChildRefIter> Elem::child_ref_range()
{
  libmesh_assert(_children);
  return {_children.get(), _children.get() + this->n_children()};
}


inline
SimpleRange<Elem::ConstChildRefIter> Elem::child_ref_range() const
{
  libmesh_assert(_children);
  return {_children.get(), _children.get() + this->n_children()};
}
#endif // LIBMESH_ENABLE_AMR




// ------------------------------------------------------------
// global Elem functions

inline
std::ostream & operator << (std::ostream & os, const Elem & e)
{
  e.print_info(os);
  return os;
}


// ------------------------------------------------------------
// Elem class member functions
inline
Elem::Elem(const unsigned int nn,
           const unsigned int ns,
           Elem * p,
           Elem ** elemlinkdata,
           Node ** nodelinkdata) :
  _nodes(nodelinkdata),
  _elemlinks(elemlinkdata),
  _sbd_id(0),
#ifdef LIBMESH_ENABLE_AMR
  _rflag(Elem::DO_NOTHING),
  _pflag(Elem::DO_NOTHING),
  _p_level(0),
#endif
  _map_type(p ? p->mapping_type() : 0),
  _map_data(p ? p->mapping_data() : 0)
{
  this->processor_id() = DofObject::invalid_processor_id;

  // If this ever legitimately fails we need to increase max_n_nodes
  libmesh_assert_less_equal(nn, max_n_nodes);

  // We currently only support refinement of elements into child
  // elements of the same type.  We can't test elem->type() here,
  // because that's virtual and we're still in the base class
  // constructor, but we can at least usually verify constency with
  // the arguments we were handed.
#ifndef NDEBUG
  if (p && !p->runtime_topology())
    {
      libmesh_assert_equal_to(nn, p->n_nodes());
      libmesh_assert_equal_to(ns, p->n_sides());
    }
#endif

  // Initialize the nodes data structure if we're given a pointer to
  // memory for it.
  if (_nodes)
    {
      for (unsigned int n=0; n<nn; n++)
        _nodes[n] = nullptr;
    }

  // Initialize the neighbors/parent data structure
  // _elemlinks = new Elem *[ns+1];

  // Initialize the elements data structure if we're given a pointer
  // to memory for it.  If we *weren't* given memory for it, e.g.
  // because a subclass like an arbitrary Polygon needs to
  // heap-allocate this memory, then that subclass will have to handle
  // this initialization too.
  if (_elemlinks)
    {
      _elemlinks[0] = p;

      for (unsigned int n=1; n<ns+1; n++)
        _elemlinks[n] = nullptr;

      // Optionally initialize data from the parent
      if (this->parent())
        {
          this->subdomain_id() = this->parent()->subdomain_id();
          this->processor_id() = this->parent()->processor_id();
          _map_type = this->parent()->_map_type;
          _map_data = this->parent()->_map_data;

#ifdef LIBMESH_ENABLE_AMR
          this->set_p_level(this->parent()->p_level());
#endif
        }
    }
}



inline
const Point & Elem::point (const unsigned int i) const
{
  libmesh_assert_less (i, this->n_nodes());
  libmesh_assert(_nodes[i]);
  libmesh_assert_not_equal_to (_nodes[i]->id(), Node::invalid_id);

  return *_nodes[i];
}



inline
Point & Elem::point (const unsigned int i)
{
  libmesh_assert_less (i, this->n_nodes());

  return *_nodes[i];
}



inline
dof_id_type Elem::node_id (const unsigned int i) const
{
  libmesh_assert_less (i, this->n_nodes());
  libmesh_assert(_nodes[i]);
  libmesh_assert_not_equal_to (_nodes[i]->id(), Node::invalid_id);

  return _nodes[i]->id();
}



inline
unsigned int Elem::local_node (const dof_id_type i) const
{
  for (auto n : make_range(this->n_nodes()))
    if (this->node_id(n) == i)
      return n;

  return libMesh::invalid_uint;
}



inline
const Node * const * Elem::get_nodes () const
{
  return _nodes;
}



inline
const Node * Elem::node_ptr (const unsigned int i) const
{
  libmesh_assert_less (i, this->n_nodes());
  libmesh_assert(_nodes[i]);

  return _nodes[i];
}



inline
Node * Elem::node_ptr (const unsigned int i)
{
  libmesh_assert_less (i, this->n_nodes());
  libmesh_assert(_nodes[i]);

  return _nodes[i];
}



inline
const Node & Elem::node_ref (const unsigned int i) const
{
  return *this->node_ptr(i);
}



inline
Node & Elem::node_ref (const unsigned int i)
{
  return *this->node_ptr(i);
}



inline
unsigned int Elem::get_node_index (const Node * node_ptr) const
{
  for (auto n : make_range(this->n_nodes()))
    if (this->_nodes[n] == node_ptr)
      return n;

  return libMesh::invalid_uint;
}



#ifdef LIBMESH_ENABLE_DEPRECATED
inline
Node * & Elem::set_node (const unsigned int i)
{
  libmesh_assert_less (i, this->n_nodes());

  libmesh_deprecated();

  return _nodes[i];
}
#endif // LIBMESH_ENABLE_DEPRECATED



inline
void Elem::set_node (const unsigned int i,
                     Node * node)
{
  libmesh_assert_less (i, this->n_nodes());

  _nodes[i] = node;
}



inline
subdomain_id_type Elem::subdomain_id () const
{
  return _sbd_id;
}



inline
subdomain_id_type & Elem::subdomain_id ()
{
  return _sbd_id;
}



inline
bool Elem::operator != (const Elem & rhs) const
{
  return !(*this == rhs);
}



inline
const Elem * Elem::neighbor_ptr (unsigned int i) const
{
  libmesh_assert_less (i, this->n_neighbors());

  return _elemlinks[i+1];
}



inline
Elem * Elem::neighbor_ptr (unsigned int i)
{
  libmesh_assert_less (i, this->n_neighbors());

  return _elemlinks[i+1];
}



inline
void Elem::set_neighbor (const unsigned int i, Elem * n)
{
  libmesh_assert_less (i, this->n_neighbors());

  _elemlinks[i+1] = n;
}



inline
bool Elem::has_neighbor (const Elem * elem) const
{
  for (auto n : this->neighbor_ptr_range())
    if (n == elem)
      return true;

  return false;
}



inline
Elem * Elem::child_neighbor (Elem * elem)
{
  for (auto n : elem->neighbor_ptr_range())
    if (n && n->parent() == this)
      return n;

  return nullptr;
}



inline
const Elem * Elem::child_neighbor (const Elem * elem) const
{
  for (auto n : elem->neighbor_ptr_range())
    if (n && n->parent() == this)
      return n;

  return nullptr;
}



inline
SimpleRange<Elem::NodeRefIter>
Elem::node_ref_range()
{
  return {_nodes, _nodes+this->n_nodes()};
}



inline
SimpleRange<Elem::ConstNodeRefIter>
Elem::node_ref_range() const
{
  return {_nodes, _nodes+this->n_nodes()};
}



inline
IntRange<unsigned short>
Elem::node_index_range() const
{
  return {0, cast_int<unsigned short>(this->n_nodes())};
}



inline
IntRange<unsigned short>
Elem::edge_index_range() const
{
  return {0, cast_int<unsigned short>(this->n_edges())};
}



inline
IntRange<unsigned short>
Elem::face_index_range() const
{
  return {0, cast_int<unsigned short>(this->n_faces())};
}



inline
IntRange<unsigned short>
Elem::side_index_range() const
{
  return {0, cast_int<unsigned short>(this->n_sides())};
}




inline
std::unique_ptr<const Elem> Elem::side_ptr (unsigned int i) const
{
  // Call the non-const version of this function, return the result as
  // a std::unique_ptr<const Elem>.
  Elem * me = const_cast<Elem *>(this);
  return me->side_ptr(i);
}



inline
void
Elem::side_ptr (std::unique_ptr<const Elem> & elem,
                const unsigned int i) const
{
  // Hand off to the non-const version of this function
  Elem * me = const_cast<Elem *>(this);
  std::unique_ptr<Elem> e {const_cast<Elem *>(elem.release())};
  me->side_ptr(e, i);
  elem = std::move(e);
}



inline
std::unique_ptr<const Elem>
Elem::build_side_ptr (const unsigned int i) const
{
  // Call the non-const version of this function, return the result as
  // a std::unique_ptr<const Elem>.
  Elem * me = const_cast<Elem *>(this);
  return me->build_side_ptr(i);
}



inline
void
Elem::build_side_ptr (std::unique_ptr<const Elem> & elem,
                      const unsigned int i) const
{
  // Hand off to the non-const version of this function
  Elem * me = const_cast<Elem *>(this);
  std::unique_ptr<Elem> e {const_cast<Elem *>(elem.release())};
  me->build_side_ptr(e, i);
  elem = std::move(e);
}



template <typename Sideclass, typename Subclass>
inline
std::unique_ptr<Elem>
Elem::simple_build_side_ptr (const unsigned int i)
{
  libmesh_assert_less (i, this->n_sides());

  std::unique_ptr<Elem> face = std::make_unique<Sideclass>();
  for (auto n : face->node_index_range())
    face->set_node(n, this->node_ptr(Subclass::side_nodes_map[i][n]));

  face->set_interior_parent(this);
  face->inherit_data_from(*this);

  return face;
}



template <typename Subclass>
inline
void
Elem::simple_build_side_ptr (std::unique_ptr<Elem> & side,
                             const unsigned int i,
                             ElemType sidetype)
{
  libmesh_assert_less (i, this->n_sides());

  if (!side.get() || side->type() != sidetype)
    {
      Subclass & real_me = cast_ref<Subclass&>(*this);
      side = real_me.Subclass::build_side_ptr(i);
    }
  else
    {
      side->set_interior_parent(this);
      side->inherit_data_from(*this);
      for (auto n : side->node_index_range())
        side->set_node(n, this->node_ptr(Subclass::side_nodes_map[i][n]));
    }
}



template <typename Subclass, typename Mapclass>
inline
void
Elem::simple_side_ptr (std::unique_ptr<Elem> & side,
                       const unsigned int i,
                       ElemType sidetype)
{
  libmesh_assert_less (i, this->n_sides());

  if (!side.get() || side->type() != sidetype)
    {
      Subclass & real_me = cast_ref<Subclass&>(*this);
      side = real_me.Subclass::side_ptr(i);
    }
  else
    {
      side->subdomain_id() = this->subdomain_id();

      for (auto n : side->node_index_range())
        side->set_node(n, this->node_ptr(Mapclass::side_nodes_map[i][n]));
    }
}



inline
std::unique_ptr<const Elem>
Elem::build_edge_ptr (const unsigned int i) const
{
  // Call the non-const version of this function, return the result as
  // a std::unique_ptr<const Elem>.
  Elem * me = const_cast<Elem *>(this);
  return me->build_edge_ptr(i);
}



inline
void
Elem::build_edge_ptr (std::unique_ptr<const Elem> & elem,
                      const unsigned int i) const
{
  // Hand off to the non-const version of this function
  Elem * me = const_cast<Elem *>(this);
  std::unique_ptr<Elem> e {const_cast<Elem *>(elem.release())};
  me->build_edge_ptr(e, i);
  elem = std::move(e);
}


template <typename Edgeclass, typename Subclass>
inline
std::unique_ptr<Elem>
Elem::simple_build_edge_ptr (const unsigned int i)
{
  libmesh_assert_less (i, this->n_edges());

  std::unique_ptr<Elem> edge = std::make_unique<Edgeclass>();

  for (auto n : edge->node_index_range())
    edge->set_node(n, this->node_ptr(Subclass::edge_nodes_map[i][n]));

  edge->set_interior_parent(this);
  edge->inherit_data_from(*this);

  return edge;
}




template <typename Subclass>
inline
void
Elem::simple_build_edge_ptr (std::unique_ptr<Elem> & edge,
                             const unsigned int i,
                             ElemType edgetype)
{
  libmesh_assert_less (i, this->n_edges());

  if (!edge.get() || edge->type() != edgetype)
    {
      Subclass & real_me = cast_ref<Subclass&>(*this);
      edge = real_me.Subclass::build_edge_ptr(i);
    }
  else
    {
      edge->inherit_data_from(*this);
      for (auto n : edge->node_index_range())
        edge->set_node(n, this->node_ptr(Subclass::edge_nodes_map[i][n]));
    }
}



inline
bool Elem::on_boundary () const
{
  // By convention, the element is on the boundary
  // if it has a nullptr neighbor.
  return this->has_neighbor(nullptr);
}



inline
unsigned int Elem::which_neighbor_am_i (const Elem * e) const
{
  libmesh_assert(e);

  const Elem * eparent = e;

  while (eparent->level() > this->level())
    {
      eparent = eparent->parent();
      libmesh_assert(eparent);
    }

  for (auto s : make_range(this->n_sides()))
    if (this->neighbor_ptr(s) == eparent)
      return s;

  return libMesh::invalid_uint;
}



inline
bool Elem::active() const
{
#ifdef LIBMESH_ENABLE_AMR
  if ((this->refinement_flag() == INACTIVE) ||
      (this->refinement_flag() == COARSEN_INACTIVE))
    return false;
  else
    return true;
#else
  return true;
#endif
}





inline
bool Elem::subactive() const
{
#ifdef LIBMESH_ENABLE_AMR
  if (this->active())
    return false;
  if (!this->has_children())
    return true;
  for (const Elem * my_ancestor = this->parent();
       my_ancestor != nullptr;
       my_ancestor = my_ancestor->parent())
    if (my_ancestor->active())
      return true;
#endif

  return false;
}



inline
bool Elem::has_children() const
{
#ifdef LIBMESH_ENABLE_AMR
  if (!_children)
    return false;
  else
    return true;
#else
  return false;
#endif
}


inline
bool Elem::has_ancestor_children() const
{
#ifdef LIBMESH_ENABLE_AMR
  if (!_children)
    return false;
  else
    for (auto & c : child_ref_range())
      if (c.has_children())
        return true;
#endif
  return false;
}



inline
bool Elem::is_ancestor_of(const Elem *
#ifdef LIBMESH_ENABLE_AMR
                          descendant
#endif
                          ) const
{
#ifdef LIBMESH_ENABLE_AMR
  const Elem * e = descendant;
  while (e)
    {
      if (this == e)
        return true;
      e = e->parent();
    }
#endif
  return false;
}



inline
const Elem * Elem::parent () const
{
  return _elemlinks[0];
}



inline
Elem * Elem::parent ()
{
  return _elemlinks[0];
}



inline
void Elem::set_parent (Elem * p)
{
  // We no longer support using parent() as interior_parent()
  libmesh_assert_equal_to(this->dim(), p ? p->dim() : this->dim());
  _elemlinks[0] = p;
}



inline
const Elem * Elem::top_parent () const
{
  const Elem * tp = this;

  // Keep getting the element's parent
  // until that parent is at level-0
  while (tp->parent() != nullptr)
    tp = tp->parent();

  libmesh_assert(tp);
  libmesh_assert_equal_to (tp->level(), 0);

  return tp;
}



inline
unsigned int Elem::level() const
{
#ifdef LIBMESH_ENABLE_AMR

  // if I don't have a parent I was
  // created directly from file
  // or by the user, so I am a
  // level-0 element
  if (this->parent() == nullptr)
    return 0;

  // if the parent and this element are of different
  // dimensionality we are at the same level as
  // the parent (e.g. we are the 2D side of a
  // 3D element)
  if (this->dim() != this->parent()->dim())
    return this->parent()->level();

  // otherwise we are at a level one
  // higher than our parent
  return (this->parent()->level() + 1);

#else

  // Without AMR all elements are
  // at level 0.
  return 0;

#endif
}



inline
unsigned int Elem::p_level() const
{
#ifdef LIBMESH_ENABLE_AMR
  return _p_level;
#else
  return 0;
#endif
}



inline
ElemMappingType Elem::mapping_type () const
{
  return static_cast<ElemMappingType>(_map_type);
}



inline
void Elem::set_mapping_type(const ElemMappingType type)
{
  _map_type = cast_int<unsigned char>(type);
}



inline
unsigned char Elem::mapping_data () const
{
  return _map_data;
}



inline
void Elem::set_mapping_data(const unsigned char data)
{
  _map_data = data;
}



#ifdef LIBMESH_ENABLE_AMR

inline
const Elem * Elem::raw_child_ptr (unsigned int i) const
{
  if (!_children)
    return nullptr;

  return _children[i];
}

inline
const Elem * Elem::child_ptr (unsigned int i) const
{
  libmesh_assert(_children);
  libmesh_assert(_children[i]);

  return _children[i];
}

inline
Elem * Elem::child_ptr (unsigned int i)
{
  libmesh_assert(_children);
  libmesh_assert(_children[i]);

  return _children[i];
}


inline
void Elem::set_child (unsigned int c, Elem * elem)
{
  libmesh_assert (this->has_children());

  _children[c] = elem;
}



inline
unsigned int Elem::which_child_am_i (const Elem * e) const
{
  libmesh_assert(e);
  libmesh_assert (this->has_children());

  unsigned int nc = this->n_children();
  for (unsigned int c=0; c != nc; c++)
    if (this->child_ptr(c) == e)
      return c;

  libmesh_error_msg("ERROR:  which_child_am_i() was called with a non-child!");

  return libMesh::invalid_uint;
}



inline
Elem::RefinementState Elem::refinement_flag () const
{
  return static_cast<RefinementState>(_rflag);
}



inline
void Elem::set_refinement_flag(RefinementState rflag)
{
  _rflag = cast_int<unsigned char>(rflag);
}



inline
Elem::RefinementState Elem::p_refinement_flag () const
{
  return static_cast<RefinementState>(_pflag);
}



inline
void Elem::set_p_refinement_flag(RefinementState pflag)
{
  if (this->p_level() == 0)
    libmesh_assert_not_equal_to
      (pflag, Elem::JUST_REFINED);

  _pflag = cast_int<unsigned char>(pflag);
}



inline
unsigned int Elem::max_descendant_p_level () const
{
  // This is undefined for subactive elements,
  // which have no active descendants
  libmesh_assert (!this->subactive());
  if (this->active())
    return this->p_level();

  unsigned int max_p_level = _p_level;
  for (auto & c : child_ref_range())
    max_p_level = std::max(max_p_level,
                           c.max_descendant_p_level());
  return max_p_level;
}



inline
void Elem::hack_p_level(unsigned int p)
{
  if (p == 0)
    libmesh_assert_not_equal_to
      (this->p_refinement_flag(), Elem::JUST_REFINED);

  _p_level = cast_int<unsigned char>(p);
}


inline
void Elem::hack_p_level_and_refinement_flag (unsigned int p,
                                             RefinementState pflag)
{
  _pflag = cast_int<unsigned char>(pflag);
  this->hack_p_level(p);
}

#endif // ifdef LIBMESH_ENABLE_AMR


inline
void Elem::orient(BoundaryInfo * boundary_info)
{
  if (this->is_flipped())
    this->flip(boundary_info);
}


inline
dof_id_type Elem::compute_key (dof_id_type n0)
{
  return n0;
}



inline
dof_id_type Elem::compute_key (dof_id_type n0,
                               dof_id_type n1)
{
  // Order the two so that n0 < n1
  if (n0 > n1) std::swap (n0, n1);

  return Utility::hashword2(n0, n1);
}



inline
dof_id_type Elem::compute_key (dof_id_type n0,
                               dof_id_type n1,
                               dof_id_type n2)
{
  std::array<dof_id_type, 3> array = {{n0, n1, n2}};
  std::sort(array.begin(), array.end());
  return Utility::hashword(array);
}



inline
dof_id_type Elem::compute_key (dof_id_type n0,
                               dof_id_type n1,
                               dof_id_type n2,
                               dof_id_type n3)
{
  std::array<dof_id_type, 4> array = {{n0, n1, n2, n3}};
  std::sort(array.begin(), array.end());
  return Utility::hashword(array);
}



inline
void Elem::inherit_data_from (const Elem & src)
{
  this->set_mapping_type(src.mapping_type());
  this->set_mapping_data(src.mapping_data());
  this->subdomain_id() = src.subdomain_id();
  this->processor_id(src.processor_id());
#ifdef LIBMESH_ENABLE_AMR
  this->set_p_level(src.p_level());
#endif
}



class Elem::SideIter
{
public:
  // Constructor with arguments.
  SideIter(const unsigned int side_number,
           Elem * parent)
    : _side(),
      _side_ptr(nullptr),
      _parent(parent),
      _side_number(side_number)
  {}


  // Empty constructor.
  SideIter()
    : _side(),
      _side_ptr(nullptr),
      _parent(nullptr),
      _side_number(libMesh::invalid_uint)
  {}


  // Copy constructor
  SideIter(const SideIter & other)
    : _side(),
      _side_ptr(nullptr),
      _parent(other._parent),
      _side_number(other._side_number)
  {}


  // op=
  SideIter & operator=(const SideIter & other)
  {
    this->_parent      = other._parent;
    this->_side_number = other._side_number;
    return *this;
  }

  // unary op*
  Elem *& operator*() const
  {
    // Set the std::unique_ptr
    this->_update_side_ptr();

    // Return a reference to _side_ptr
    return this->_side_ptr;
  }

  // op++
  SideIter & operator++()
  {
    ++_side_number;
    return *this;
  }

  // op==  Two side iterators are equal if they have
  // the same side number and the same parent element.
  bool operator == (const SideIter & other) const
  {
    return (this->_side_number == other._side_number &&
            this->_parent      == other._parent);
  }


  // Consults the parent Elem to determine if the side
  // is a boundary side.  Note: currently side N is a
  // boundary side if neighbor N is nullptr.  Be careful,
  // this could possibly change in the future?
  bool side_on_boundary() const
  {
    return this->_parent->neighbor_ptr(_side_number) == nullptr;
  }

private:
  // Update the _side pointer by building the correct side.
  // This has to be called before dereferencing.
  void _update_side_ptr() const
  {
    // Construct new side, store in std::unique_ptr
    this->_side = this->_parent->build_side_ptr(this->_side_number);

    // Also set our internal naked pointer.  Memory is still owned
    // by the std::unique_ptr.
    this->_side_ptr = _side.get();
  }

  // std::unique_ptr to the actual side, handles memory management for
  // the sides which are created during the course of iteration.
  mutable std::unique_ptr<Elem> _side;

  // Raw pointer needed to facilitate passing back to the user a
  // reference to a non-temporary raw pointer in order to conform to
  // the variant_filter_iterator interface.  It points to the same
  // thing the std::unique_ptr "_side" above holds.  What happens if the user
  // calls delete on the pointer passed back?  Well, this is an issue
  // which is not addressed by the iterators in libMesh.  Basically it
  // is a bad idea to ever call delete on an iterator from the library.
  mutable Elem * _side_ptr;

  // Pointer to the parent Elem class which generated this iterator
  Elem * _parent;

  // A counter variable which keeps track of the side number
  unsigned int _side_number;
};






// Private implementation functions in the Elem class for the side iterators.
// They have to come after the definition of the SideIter class.
inline
Elem::SideIter Elem::_first_side()
{
  return SideIter(0, this);
}



inline
Elem::SideIter Elem::_last_side()
{
  return SideIter(this->n_neighbors(), this);
}




struct
Elem::side_iterator : variant_filter_iterator<Elem::Predicate, Elem *>
{
  // Templated forwarding ctor -- forwards to appropriate variant_filter_iterator ctor
  template <typename PredType, typename IterType>
  side_iterator (const IterType & d,
                 const IterType & e,
                 const PredType & p ) :
    variant_filter_iterator<Elem::Predicate, Elem *>(d,e,p) {}
};



inline
SimpleRange<Elem::NeighborPtrIter> Elem::neighbor_ptr_range()
{
  return {_elemlinks+1, _elemlinks + 1 + this->n_neighbors()};
}


inline
SimpleRange<Elem::ConstNeighborPtrIter> Elem::neighbor_ptr_range() const
{
  return {_elemlinks+1, _elemlinks + 1 + this->n_neighbors()};
}

} // namespace libMesh


// Helper function for default caches in Elem subclasses

#define LIBMESH_ENABLE_TOPOLOGY_CACHES                                  \
  virtual                                                               \
  std::vector<std::vector<std::vector<std::vector<std::pair<unsigned char, unsigned char>>>>> & \
  _get_bracketing_node_cache() const override                   \
  {                                                                     \
    static std::vector<std::vector<std::vector<std::vector<std::pair<unsigned char, unsigned char>>>>> c; \
    return c;                                                           \
  }                                                                     \
                                                                        \
  virtual                                                               \
  std::vector<std::vector<std::vector<signed char>>> &                  \
  _get_parent_indices_cache() const override                    \
  {                                                                     \
    static std::vector<std::vector<std::vector<signed char>>> c;        \
    return c;                                                           \
  }






#endif // LIBMESH_ELEM_H