elem.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2019 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
 
#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;
 
  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);
 
  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);
 
  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_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;
 
  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;
 
  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;
 
  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 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;
 
  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 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);
 
  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);
 
  template <typename Subclass, typename Mapclass>
  void simple_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
 
  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, _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
  _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);
 
  // 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();
      _map_type = this->parent()->_map_type;
      _map_data = this->parent()->_map_data;
    }
 
#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();
}
 
 
 
inline
unsigned int Elem::local_node (const dof_id_type i) const
{
  for (auto n : IntRange<unsigned int>(0, 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 : IntRange<unsigned int>(0, this->n_nodes()))
    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];
}
 
 
 
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());
}
 
 
 
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());
}
 
 
 
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, 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]);
    }
}
 
 
 
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);
  const Elem * e = const_cast<const Elem *>(me->build_edge_ptr(i).release());
  return std::unique_ptr<const Elem>(e);
}
 
 
 
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 : IntRange<unsigned int>(0, 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 == 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)
{
  // 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);
}
 
 
 
#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