elem.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 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/dof_object.h"
#include "libmesh/id_types.h"
#include "libmesh/reference_counted_object.h"
#include "libmesh/node.h"
#include "libmesh/enum_elem_type.h"
#include "libmesh/enum_elem_quality.h"
#include "libmesh/enum_order.h"
#include "libmesh/enum_io_package.h"
#include "libmesh/auto_ptr.h"
#include "libmesh/multi_predicates.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>

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:

  virtual ~Elem();

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

  Point & point (const unsigned int i);

  dof_id_type node (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;

  Node* get_node (const unsigned int i) const;

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

  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;

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

  Elem* neighbor (const unsigned int i) 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,
                                 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;

  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;

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

  //   /**
  //    * @returns the VTK element type of the sc-th sub-element.
  //    */
  //   virtual unsigned int vtk_element_type (const unsigned int sc) const = 0;

  virtual ElemType type () const = 0;

  virtual unsigned int 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 type_to_n_sides_map[INVALID_ELEM];

  virtual unsigned int n_sides () const = 0;

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

  virtual unsigned int n_vertices () const = 0;

  virtual unsigned int n_edges () const = 0;

  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 bool is_edge(const unsigned int i) const = 0;

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

  /*
   * @returns true iff the specified (local) node number is on the
   * specified side
   */
  virtual bool is_node_on_side(const unsigned int n,
                               const unsigned int s) const = 0;

  /*
   * @returns true iff the specified (local) node number is on the
   * specified edge
   */
  virtual bool is_node_on_edge(const unsigned int n,
                               const unsigned int e) const = 0;

  /*
   * @returns true iff the specified edge is on the specified side
   */
  virtual bool is_edge_on_side(const unsigned int e,
                               const unsigned int s) const = 0;

  /*
   * @returns the side number opposite to \p s (for a tensor product
   * element), or throws an error otherwise.
   */
  virtual unsigned int opposite_side(const unsigned int s) const;

  /*
   * @returns the local node number for the node opposite to node n
   * on side \p opposite_side(s) (for a tensor product element), or
   * throws an error otherwise.
   */
  virtual unsigned int opposite_node(const unsigned int n,
                                     const unsigned int s) const;

  //   /**
  //    * @returns the number of children this element has that
  //    * share side \p s
  //    */
  //   virtual unsigned int n_children_per_side (const unsigned int) const = 0;

  virtual unsigned int n_sub_elem () const = 0;

  virtual UniquePtr<Elem> side (const unsigned int i) const = 0;

  virtual UniquePtr<Elem> build_side (const unsigned int i,
                                      bool proxy=true) const = 0;

  virtual UniquePtr<Elem> build_edge (const unsigned int i) const = 0;

  virtual Order default_order () const = 0;

  virtual Point centroid () const;

  virtual Real hmin () const;

  virtual Real hmax () const;

  virtual Real volume () 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;

  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;

#ifdef LIBMESH_ENABLE_AMR

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

  Elem* child (const unsigned int i) 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_side(const unsigned int c,
                                const unsigned int s) const = 0;

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

  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;
  //typedef variant_filter_iterator<Elem*, Predicate> side_iterator;

  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 Point origin () const { libmesh_not_implemented(); return Point(); }

#endif




  static UniquePtr<Elem> build (const ElemType type,
                                Elem* p=NULL);

#ifdef LIBMESH_ENABLE_AMR

  virtual float embedding_matrix (const unsigned int i,
                                  const unsigned int j,
                                  const unsigned int k) const = 0;

#endif


protected:

  //-------------------------------------------------------
  // These methods compute has keys from the specified
  // global node numbers
  //
  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);
  //-------------------------------------------------------



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;
  //RefinementState _rflag;

  unsigned char _pflag;
  //RefinementState _pflag;

  unsigned char _p_level;
#endif
};

// ------------------------------------------------------------
// 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(NULL),
#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;

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

  // 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] = NULL;
    }

  // Optionally initialize data from the parent
  if (this->parent() != NULL)
    {
      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 != NULL)
  //   delete [] _nodes;
  // _nodes = NULL;

  // delete [] _elemlinks;

#ifdef LIBMESH_ENABLE_AMR

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

#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 (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 (unsigned int n=0; n != this->n_nodes(); ++n)
    if (this->node(n) == i)
      return n;

  return libMesh::invalid_uint;
}



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



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

  return _nodes[i];
}



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
Elem* Elem::neighbor (const unsigned int i) const
{
  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 (unsigned int n=0; n<this->n_neighbors(); n++)
    if (this->neighbor(n) == elem)
      return true;

  return false;
}



inline
Elem* Elem::child_neighbor (Elem* elem) const
{
  for (unsigned int n=0; n<elem->n_neighbors(); n++)
    if (elem->neighbor(n) &&
        elem->neighbor(n)->parent() == this)
      return elem->neighbor(n);

  return NULL;
}



inline
const Elem* Elem::child_neighbor (const Elem* elem) const
{
  for (unsigned int n=0; n<elem->n_neighbors(); n++)
    if (elem->neighbor(n) &&
        elem->neighbor(n)->parent() == this)
      return elem->neighbor(n);

  return NULL;
}



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



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; s<this->n_neighbors(); s++)
    if (this->neighbor(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
      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 != NULL;
       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 == NULL)
    return false;
  else
    return true;
#else
  return false;
#endif
}


inline
bool Elem::has_ancestor_children() const
{
#ifdef LIBMESH_ENABLE_AMR
  if (_children == NULL)
    return false;
  else
    for (unsigned int c=0; c != this->n_children(); c++)
      if (this->child(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() != NULL)
    tp = tp->parent();

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

  return tp;
}



inline
const Elem* Elem::interior_parent () const
{
  // interior parents make no sense for full-dimensional elements.
  libmesh_assert_less (this->dim(), LIBMESH_DIM);

  // // and they [USED TO BE] only good for level-0 elements
  // if (this->level() != 0)
  // return this->parent()->interior_parent();

  // We store the interior_parent pointer after both the parent
  // neighbor and neighbor pointers
  Elem *interior_p = _elemlinks[1+this->n_sides()];

  // If we have an interior_parent, it had better be the
  // one-higher-dimensional interior element we are looking for.
  libmesh_assert (!interior_p ||
                  interior_p->dim() == (this->dim()+1));

  return interior_p;
}



inline
void Elem::set_interior_parent (Elem *p)
{
  // interior parents make no sense for full-dimensional elements.
  libmesh_assert_less (this->dim(), LIBMESH_DIM);

  // this had better be a one-higher-dimensional interior element
  libmesh_assert (!p ||
                  p->dim() == (this->dim()+1));

  _elemlinks[1+this->n_sides()] = p;
}



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() == NULL)
    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
Elem* Elem::child (const unsigned int i) const
{
  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());

  for (unsigned int c=0; c<this->n_children(); c++)
    if (this->child(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)
{
  _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 (unsigned int c=0; c != this->n_children(); c++)
    max_p_level = std::max(max_p_level,
                           this->child(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() != NULL)
    {
      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);
        }
      // 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 (unsigned int c=0; c != this->parent()->n_children(); c++)
            parent_p_level = std::min(parent_p_level,
                                      this->parent()->child(c)->p_level());

          if (parent_p_level != this->parent()->p_level())
            this->parent()->set_p_level(parent_p_level);

          return;
        }
    }

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



inline
void Elem::hack_p_level(unsigned int p)
{
  _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(NULL),
      _parent(parent),
      _side_number(side_number)
  {}


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


  // Copy constructor
  SideIter(const SideIter& other)
    : _side(),
      _side_ptr(NULL),
      _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 UniquePtr
    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 nieghbor N is NULL.  Be careful,
  // this could possibly change in the future?
  bool side_on_boundary() const
  {
    return this->_parent->neighbor(_side_number) == NULL;
  }

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 UniquePtr
    this->_side = this->_parent->build_side(this->_side_number);

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

  // UniquePtr to the actual side, handles memory management for
  // the sides which are created during the course of iteration.
  mutable UniquePtr<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 UniquePtr "_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) {}
};


} // namespace libMesh


#endif // LIBMESH_ELEM_H