elem.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2018 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/auto_ptr.h" // deprecated
#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

#ifdef LIBMESH_FORWARD_DECLARE_ENUMS
namespace libMesh
{
enum ElemQuality : int;
enum IOPackage : int;
enum Order : int;
}
#else
#include "libmesh/enum_elem_quality.h"
#include "libmesh/enum_io_package.h"
#include "libmesh/enum_order.h"
#endif

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

namespace libMesh
{

// Forward declarations
class MeshBase;
class MeshRefinement;
class Elem;
#ifdef LIBMESH_ENABLE_PERIODIC
class PeriodicBoundaries;
class PointLocatorBase;
#endif

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();

  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;

#ifdef LIBMESH_ENABLE_DEPRECATED
  dof_id_type node (const unsigned int i) const;
#endif

  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
  Node * get_node (const unsigned int i) const;
#endif

  virtual Node * & set_node (const unsigned int i);

  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 const subdomain_id_type invalid_subdomain_id;

  const Elem * reference_elem () const;

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

  virtual dof_id_type key () const;

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

  const Elem * neighbor_ptr (unsigned int i) const;

  Elem * neighbor_ptr (unsigned int i);

#ifdef LIBMESH_ENABLE_DEPRECATED
  Elem * neighbor (const unsigned int i) const;
#endif

  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 which_node_am_i(unsigned int side,
                                       unsigned int side_node) const = 0;

  bool contains_vertex_of(const Elem * e) 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_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 remove_links_to_me ();

  void make_links_to_me_remote ();

  void make_links_to_me_local (unsigned int n);

  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;

  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;

  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;

  virtual unsigned int n_children () const = 0;

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

  virtual unsigned int 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;

  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 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;


#ifdef LIBMESH_ENABLE_DEPRECATED
  std::unique_ptr<Elem> side (const unsigned int i) const;
#endif

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

  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;

#ifdef LIBMESH_ENABLE_DEPRECATED
  std::unique_ptr<Elem> build_side (const unsigned int i, bool proxy=true) const;
#endif

  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;

#ifdef LIBMESH_ENABLE_DEPRECATED
  std::unique_ptr<Elem> build_edge (const unsigned int i) const;
#endif

  virtual Order default_order () const = 0;

  virtual Point centroid () const;

  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 close_to_point(const Point & p, Real tol) const;

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 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);


  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;

#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);

#ifdef LIBMESH_ENABLE_DEPRECATED
  Elem * child (const unsigned int i) const;
#endif

  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,
                    const bool reset=true) const;

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

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

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

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

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

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

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

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

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

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

  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);

  virtual void refine (MeshRefinement & mesh_refinement);

  void coarsen ();

  void contract ();

#endif

#ifdef DEBUG

  void libmesh_assert_valid_neighbors() const;

  void libmesh_assert_valid_node_pointers() const;
#endif // DEBUG

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(); }

#endif




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

#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 float 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 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);

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

#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

  Elem ** _children;
#endif

  subdomain_id_type _sbd_id;

#ifdef LIBMESH_ENABLE_AMR

  unsigned char _rflag;

  unsigned char _pflag;

  unsigned char _p_level;
#endif
};



// ------------------------------------------------------------
// 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, _children + this->n_children()};
}


inline
SimpleRange<Elem::ConstChildRefIter> Elem::child_ref_range() const
{
  libmesh_assert(_children);
  return {_children, _children + 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),
#ifdef LIBMESH_ENABLE_AMR
  _children(nullptr),
#endif
  _sbd_id(0)
#ifdef LIBMESH_ENABLE_AMR
  ,
  _rflag(Elem::DO_NOTHING),
  _pflag(Elem::DO_NOTHING),
  _p_level(0)
#endif
{
  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);

  // Initialize the nodes data structure
  if (_nodes)
    {
      for (unsigned int n=0; n<nn; n++)
        _nodes[n] = nullptr;
    }

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

  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() != nullptr)
    {
      this->subdomain_id() = this->parent()->subdomain_id();
      this->processor_id() = this->parent()->processor_id();
    }

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



inline
Elem::~Elem()
{
  // Deleting my parent/neighbor/nodes storage isn't necessary since it's
  // handled by the subclass

  // if (_nodes != nullptr)
  //   delete [] _nodes;
  // _nodes = nullptr;

  // delete [] _elemlinks;

#ifdef LIBMESH_ENABLE_AMR

  // Delete my children's storage
  if (_children != nullptr)
    delete [] _children;
  _children = nullptr;

#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();
}



#ifdef LIBMESH_ENABLE_DEPRECATED
inline
dof_id_type Elem::node (const unsigned int i) const
{
  libmesh_deprecated();
  return this->node_id(i);
}
#endif



inline
unsigned int Elem::local_node (const dof_id_type i) const
{
  for (unsigned int n=0; n != this->n_nodes(); ++n)
    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);
}



#ifdef LIBMESH_ENABLE_DEPRECATED
inline
Node * Elem::get_node (const unsigned int i) const
{
  // This const function has incorrectly returned a non-const pointer
  // for years.  Now that it is reimplemented in terms of the new
  // interface which does return a const pointer, we need to use a
  // const_cast to mimic the old (incorrect) behavior.  This function
  // is now deprecated and eventually will be removed entirely,
  // obviating the need for this ugly cast.
  libmesh_deprecated();
  return const_cast<Node *>(this->node_ptr(i));
}
#endif



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

  return libMesh::invalid_uint;
}



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

  return _nodes[i];
}



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



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



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];
}



#ifdef LIBMESH_ENABLE_DEPRECATED
inline
Elem * Elem::neighbor (const unsigned int i) const
{
  // Support the deprecated interface by calling the new,
  // const-correct interface and casting the result to an Elem *.
  libmesh_deprecated();
  return const_cast<Elem *>(this->neighbor_ptr(i));
}
#endif



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::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);
  const Elem * s = const_cast<const Elem *>(me->side_ptr(i).release());
  return std::unique_ptr<const Elem>(s);
}



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.reset(e.release());
}



#ifdef LIBMESH_ENABLE_DEPRECATED
inline
std::unique_ptr<Elem> Elem::side (const unsigned int i) const
{
  // Call the const version of side_ptr(), and const_cast the result.
  libmesh_deprecated();
  Elem * s = const_cast<Elem *>(this->side_ptr(i).release());
  return std::unique_ptr<Elem>(s);
}
#endif



inline
std::unique_ptr<const Elem>
Elem::build_side_ptr (const unsigned int i, bool proxy) 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);
  const Elem * s = const_cast<const Elem *>(me->build_side_ptr(i, proxy).release());
  return std::unique_ptr<const Elem>(s);
}



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.reset(e.release());
}



#ifdef LIBMESH_ENABLE_DEPRECATED
inline
std::unique_ptr<Elem>
Elem::build_side (const unsigned int i, bool proxy) const
{
  // Call the const version of build_side_ptr(), and const_cast the result.
  libmesh_deprecated();
  Elem * s = const_cast<Elem *>(this->build_side_ptr(i, proxy).release());
  return std::unique_ptr<Elem>(s);
}
#endif



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)
    side = this->build_side_ptr(i, false);
  else
    {
      side->subdomain_id() = this->subdomain_id();

      for (auto n : side->node_index_range())
        side->set_node(n) = this->node_ptr(Subclass::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);
  const Elem * e = const_cast<const Elem *>(me->build_edge_ptr(i).release());
  return std::unique_ptr<const Elem>(e);
}



#ifdef LIBMESH_ENABLE_DEPRECATED
inline
std::unique_ptr<Elem>
Elem::build_edge (const unsigned int i) const
{
  // Call the const version of build_edge_ptr(), and const_cast the result.
  libmesh_deprecated();
  Elem * e = const_cast<Elem *>(this->build_edge_ptr(i).release());
  return std::unique_ptr<Elem>(e);
}
#endif



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 (unsigned int s=0, n_s = this->n_sides(); s != n_s; ++s)
    if (this->neighbor_ptr(s) == eparent)
      return s;

  return libMesh::invalid_uint;
}



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

  const unsigned int ns = this->n_sides();
  const unsigned int nn = this->n_nodes();

  const unsigned int en = e->n_nodes();

  // e might be on any side until proven otherwise
  std::vector<bool> might_be_side(ns, true);

  for (unsigned int i=0; i != en; ++i)
    {
      Point side_point = e->point(i);
      unsigned int local_node_id = libMesh::invalid_uint;

      // Look for a node of this that's contiguous with node i of
      // e. Note that the exact floating point comparison of Point
      // positions is intentional, see the class documentation for
      // this function.
      for (unsigned int j=0; j != nn; ++j)
        if (this->point(j) == side_point)
          local_node_id = j;

      // If a node of e isn't contiguous with some node of this, then
      // e isn't a side of this.
      if (local_node_id == libMesh::invalid_uint)
        return libMesh::invalid_uint;

      // If a node of e isn't contiguous with some node on side s of
      // this, then e isn't on side s.
      for (unsigned int s=0; s != ns; ++s)
        if (!this->is_node_on_side(local_node_id, s))
          might_be_side[s] = false;
    }

  for (unsigned int s=0; s != ns; ++s)
    if (might_be_side[s])
      {
#ifdef DEBUG
        for (unsigned int s2=s+1; s2 < ns; ++s2)
          libmesh_assert (!might_be_side[s2]);
#endif
        return s;
      }

  // Didn't find any matching side
  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 == nullptr)
    return false;
  else
    return true;
#else
  return false;
#endif
}


inline
bool Elem::has_ancestor_children() const
{
#ifdef LIBMESH_ENABLE_AMR
  if (_children == nullptr)
    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)
{
  _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
}



#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];
}


#ifdef LIBMESH_ENABLE_DEPRECATED
inline
Elem * Elem::child (const unsigned int i) const
{
  // Support the deprecated interface by calling the new,
  // const-correct interface and casting the result to an Elem *.
  libmesh_deprecated();
  return const_cast<Elem *>(this->child_ptr(i));
}
#endif



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<RefinementState>(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::set_p_level(unsigned int p)
{
  // Maintain the parent's p level as the minimum of it's children
  if (this->parent() != nullptr)
    {
      unsigned int parent_p_level = this->parent()->p_level();

      // If our new p level is less than our parents, our parents drops
      if (parent_p_level > p)
        {
          this->parent()->set_p_level(p);

          // And we should keep track of the drop, in case we need to
          // do a projection later.
          this->parent()->set_p_refinement_flag(Elem::JUST_COARSENED);
        }
      // If we are the lowest p level and it increases, so might
      // our parent's, but we have to check every other child to see
      else if (parent_p_level == _p_level && _p_level < p)
        {
          _p_level = cast_int<unsigned char>(p);
          parent_p_level = cast_int<unsigned char>(p);
          for (auto & c : this->parent()->child_ref_range())
            parent_p_level = std::min(parent_p_level,
                                      c.p_level());

          // When its children all have a higher p level, the parent's
          // should rise
          if (parent_p_level > this->parent()->p_level())
            {
              this->parent()->set_p_level(parent_p_level);

              // And we should keep track of the rise, in case we need to
              // do a projection later.
              this->parent()->set_p_refinement_flag(Elem::JUST_REFINED);
            }

          return;
        }
    }

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



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);
}



#endif // ifdef LIBMESH_ENABLE_AMR


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)
{
  // Order the numbers such that n0 < n1 < n2.
  // We'll do it in 3 steps like this:
  //
  //     n0         n1                n2
  //     min(n0,n1) max(n0,n1)        n2
  //     min(n0,n1) min(n2,max(n0,n1) max(n2,max(n0,n1)
  //           |\   /|                  |
  //           | \ / |                  |
  //           |  /  |                  |
  //           | /  \|                  |
  //  gb min= min   max              gb max

  // Step 1
  if (n0 > n1) std::swap (n0, n1);

  // Step 2
  if (n1 > n2) std::swap (n1, n2);

  // Step 3
  if (n0 > n1) std::swap (n0, n1);

  libmesh_assert ((n0 < n1) && (n1 < n2));

  dof_id_type array[3] = {n0, n1, n2};
  return Utility::hashword(array, 3);
}



inline
dof_id_type Elem::compute_key (dof_id_type n0,
                               dof_id_type n1,
                               dof_id_type n2,
                               dof_id_type n3)
{
  // Sort first
  // Step 1
  if (n0 > n1) std::swap (n0, n1);

  // Step 2
  if (n2 > n3) std::swap (n2, n3);

  // Step 3
  if (n0 > n2) std::swap (n0, n2);

  // Step 4
  if (n1 > n3) std::swap (n1, n3);

  // Finally sort step 5
  if (n1 > n2) std::swap (n1, n2);

  libmesh_assert ((n0 < n1) && (n1 < n2) && (n2 < n3));

  dof_id_type array[4] = {n0, n1, n2, n3};
  return Utility::hashword(array, 4);
}



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