elem.h

一个用来实现偏微分方程中网格的计算库

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



inline
void Elem::set_p_refinement_flag(RefinementState pflag)
{
#ifdef DEBUG
  if (pflag != static_cast<RefinementState>(static_cast<unsigned char>(pflag)))
    {
      std::cerr << "ERROR: unsigned char too small to hold Elem::_pflag!"
		<< std::endl
		<< "Recompile with Elem:_*flag set to something bigger!"
		<< std::endl;
      libmesh_error();
    }
#endif

  _pflag = 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)
{
#ifdef DEBUG
  if (p != static_cast<unsigned int>(static_cast<unsigned char>(p)))
    {
      std::cerr << "ERROR: unsigned char too small to hold Elem::_p_level!"
		<< std::endl
		<< "Recompile with Elem:_p_level set to something bigger!"
		<< std::endl;
      libmesh_error();
    }
#endif

  // 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 = p;
	  parent_p_level = 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 = p;
}



inline
void Elem::hack_p_level(unsigned int p)
{
  _p_level = p;
}



#endif /* ifdef ENABLE_AMR */


inline
unsigned int Elem::compute_key (unsigned int n0)
{
  return n0;
}



inline
unsigned int Elem::compute_key (unsigned int n0,
				unsigned int n1)
{
  // big prime number
  const unsigned int bp = 65449;
  
  // Order the two so that n0 < n1
  if (n0 > n1) std::swap (n0, n1);

  return (n0%bp + (n1<<5)%bp);  
}



inline
unsigned int Elem::compute_key (unsigned int n0,
				unsigned int n1,
				unsigned int n2)
{
  // big prime number
  const unsigned int bp = 65449;

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

  
  return (n0%bp + (n1<<5)%bp + (n2<<10)%bp);
}



inline
unsigned int Elem::compute_key (unsigned int n0,
				unsigned int n1,
				unsigned int n2,
				unsigned int n3)
{
  // big prime number
  const unsigned int bp = 65449;

  // 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 step 5
  if (n1 > n2) std::swap (n1, n2);

  libmesh_assert ((n0 < n1) && (n1 < n2) && (n2 < n3));
  
  return (n0%bp + (n1<<5)%bp + (n2<<10)%bp + (n3<<15)%bp);
}



//-----------------------------------------------------------------
/**
 * Convenient way to communicate elements.  This class
 * packes up an element so that it can easily be communicated through
 * an MPI array.
 *
 * \author Benjamin S. Kirk
 * \date 2008
 */
class Elem::PackedElem
{
private:
    
  /**
   * Iterator pointing to the beginning of this packed element's index buffer.
   */
  const std::vector<int>::const_iterator _buf_begin;

public:

  /**
   * Constructor.  Takes an input iterator pointing to the 
   * beginning of the connectivity block for this element.
   */
  PackedElem (const std::vector<int>::const_iterator _buf_in) :
    _buf_begin(_buf_in)
  {}
    
  /**
   * An \p Elem can be packed into an integer array which is 
   * \p header_size + elem->n_nodes() in length.
   */
  static const unsigned int header_size; /* = 10 with AMR, 4 without */
    
  /**
   * For each element it is of the form
   * [ level p_level r_flag p_r_flag etype processor_id subdomain_id 
   *  self_ID parent_ID which_child node_0 node_1 ... node_n]
   * We cannot use unsigned int because parent_ID can be negative
   */
  static void pack (std::vector<int> &conn, const Elem* elem);    

  /**
   * Unpacks this packed element.  Returns a pointer to the new element.
   * Takes a pointer to the parent, which is required unless this packed
   * element is at level 0.
   */
  Elem * unpack (MeshBase &mesh, Elem *parent = NULL) const;

  /**
   * \p return the level of this packed element.
   */
  unsigned int level () const
  {
    return static_cast<unsigned int>(*_buf_begin);
  }

  /**
   * \p return the p-level of this packed element.
   */
  unsigned int p_level () const
  {
    return static_cast<unsigned int>(*(_buf_begin+1));
  }

#ifdef ENABLE_AMR
  /**
   * \p return the refinement state of this packed element.
   */
  Elem::RefinementState refinement_flag () const
  {
    return static_cast<Elem::RefinementState>(*(_buf_begin+2));
  }

  /**
   * \p return the p-refinement state of this packed element.
   */
  Elem::RefinementState p_refinement_flag () const
  {
    return static_cast<Elem::RefinementState>(*(_buf_begin+3));
  }
#endif // ENABLE_AMR

  /**
   * \p return the element type of this packed element.
   */
  ElemType type () const
  {
    return static_cast<ElemType>(*(_buf_begin+4));
  }

  /**
   * \p return the processor id of this packed element.
   */
  unsigned int processor_id () const
  {
    return static_cast<unsigned int>(*(_buf_begin+5));
  }

  /**
   * \p return the subdomain id of this packed element.
   */
  unsigned int subdomain_id () const
  {
    return static_cast<unsigned int>(*(_buf_begin+6));
  }

  /**
   * \p return the id of this packed element.
   */
  unsigned int id () const
  {
    return static_cast<unsigned int>(*(_buf_begin+7));
  }

  /**
   * \p return the parent id of this packed element.
   */
  int parent_id () const
  {
    return *(_buf_begin+8);
  }

  /**
   * \p return which child this packed element is.
   */
  unsigned int which_child_am_i () const
  {
    return static_cast<unsigned int>(*(_buf_begin+9));
  }
    
  /**
   * \p return the number of nodes in this packed element
   */
  unsigned int n_nodes () const
  {
    return Elem::type_to_n_nodes_map[this->type()];
  }

  /**
   * \p return the global index of the packed element's nth node.
   */
  unsigned int node (const unsigned int n) const
  {
    return static_cast<unsigned int>(*(_buf_begin+10+n));
  }    
}; // end class PackedElem


/**
 * The definition of the protected nested SideIter class.
 */
class Elem::SideIter
{
public:
  // Constructor with arguments.
  SideIter(const unsigned int side_number,
	   Elem* parent)
    : _side_number(side_number),
      _side_ptr(NULL),
      _parent(parent)
  {}

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


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


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

  // unary op*
  Elem*& operator*() const
  {
    // Set the AutoPtr
    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 AutoPtr
    this->_side = this->_parent->build_side(this->_side_number);

    // Also set our internal naked pointer.  Memory is still owned
    // by the AutoPtr.
    this->_side_ptr = _side.get();
  }
    
  // A counter variable which keeps track of the side number
  unsigned int _side_number;
    
  // AutoPtr to the actual side, handles memory management for
  // the sides which are created during the course of iteration.
  mutable AutoPtr<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 AutoPtr "_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;

};






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



/**
 * The definition of the struct used for iterating over sides.
 */
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) {}
};




#endif // end #ifndef __elem_h__