dof_map.h

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

   * Constrains the element vector.
   */
  void constrain_element_vector (DenseVector<Number>&       rhs,
				 std::vector<unsigned int>& dofs,
				 bool asymmetric_constraint_rows = true) const;
  
  /**
   * Constrains the element matrix and vector.  This method requires
   * the element matrix to be square, in which case the elem_dofs
   * correspond to the global DOF indices of both the rows and
   * columns of the element matrix.  For this case the rows
   * and columns of the matrix necessarily correspond to variables
   * of the same approximation order.
   */
  void constrain_element_matrix_and_vector (DenseMatrix<Number>& matrix,
					    DenseVector<Number>& rhs,
					    std::vector<unsigned int>& elem_dofs,
					    bool asymmetric_constraint_rows = true) const;
  
  /**
   * Constrains the numeric vector \p v, which represents a solution defined on
   * the mesh.  This may need to be used after a linear solve, if your linear
   * solver's solutions do not satisfy your DoF constraints to a tight enough
   * tolerance.
   *
   * If \p v == NULL, the system solution vector is constrained
   */
  void enforce_constraints_exactly (const System &system,
				    NumericVector<Number> *v = NULL) const;
  
  /**
   * Prints the \p _dof_constraints data structure.
   */
  void print_dof_constraints(std::ostream& os=std::cout) const;

  /**
   * Tests the constrained degrees of freedom on the numeric vector \p v, which
   * represents a solution defined on the mesh, returning a pair whose first
   * entry is the maximum absolute error on a constrained DoF and whose second
   * entry is the maximum relative error.  Useful for debugging purposes.
   *
   * If \p v == NULL, the system solution vector is tested.
   */
  std::pair<Real, Real> max_constraint_error(const System &system,
				    NumericVector<Number> *v = NULL) const;

#endif // ENABLE_AMR || ENABLE_PERIODIC


#ifdef ENABLE_PERIODIC

  //--------------------------------------------------------------------
  // PeriodicBoundary-specific methods

  /**
   * Adds a copy of the specified periodic boundary to the system.
   */
  void add_periodic_boundary (const PeriodicBoundary& periodic_boundary);

  /**
   * @returns true if the boundary given by \p boundaryid is periodic,
   * false otherwise
   */
  bool is_periodic_boundary (const unsigned int boundaryid) const;

#endif // ENABLE_PERIODIC


#ifdef ENABLE_AMR

  //--------------------------------------------------------------------
  // AMR-specific methods

  /**
   * After a mesh is refined and repartitioned it is possible that the
   * \p _send_list will need to be augmented.  This is the case when an
   * element is refined and its children end up on different processors
   * than the parent.  These children will need values from the parent
   * when projecting the solution onto the refined mesh, hence the parent's
   * DOF indices need to be included in the \p _send_list.
   */
  void augment_send_list_for_projection(const MeshBase &);

  /**
   * Fills the vector di with the global degree of freedom indices
   * for the element using the \p DofMap::old_dof_object.
   * If no variable number is specified then all
   * variables are returned.
   */
  void old_dof_indices (const Elem* const elem,
			std::vector<unsigned int>& di,
			const unsigned int vn = libMesh::invalid_uint) const;
  /**
   * @returns the total number of degrees of freedom on old_dof_objects
   * 
   */
  unsigned int n_old_dofs() const { return _n_old_dfs; }

  /**
   * Constrains degrees of freedom on side \p s of element \p elem which
   * correspond to variable number \p var and to p refinement levels 
   * above \p p.
   */
  void constrain_p_dofs (unsigned int var,
			 const Elem *elem,
			 unsigned int s,
			 unsigned int p);
  
#endif // ENABLE_AMR

  /**
   * Reinitialize the underlying data strucures conformal to the current mesh.
   */
  void reinit (MeshBase& mesh);
  
  /**
   * Free all memory associated with the object, but keep the mesh pointer.
   */
  void clear ();
  
  /**
   * Degree of freedom coupling.  If left empty each DOF
   * couples to all others.  Can be used to reduce memory
   * requirements for sparse matrices.  DOF 0 might only
   * couple to itself, in which case \p dof_coupling(0,0)
   * should be 1 and \p dof_coupling(0,j) = 0 for j not equal
   * to 0.
   *
   * This variable is named as though it were class private,
   * but it is in the public interface.  Also there are no
   * public methods for accessing it...  This typically means
   * you should only use it if you know what you are doing.
   */
  CouplingMatrix* _dof_coupling;

  
private:
  
  /**
   * @returns the number of the system we are responsible for.
   */
  unsigned int sys_number() const;

  /**
   * @invalidates all active DofObject dofs for this system
   */
  void invalidate_dofs(MeshBase& mesh) const;

  /**
   * An adapter function that returns Node pointers by index
   */
  DofObject* node_ptr(MeshBase& mesh, unsigned int i) const;

  /**
   * An adapter function that returns Elem pointers by index
   */
  DofObject* elem_ptr(MeshBase& mesh, unsigned int i) const;

  /**
   * A member function type like node_ptr or elem_ptr
   */
  typedef DofObject* (DofMap::*dofobject_accessor)
    (MeshBase& mesh, unsigned int i) const;
  
  /**
   * Helper function for distributing dofs in parallel
   */
  template<typename iterator_type>
  void set_nonlocal_dof_objects(iterator_type objects_begin,
			        iterator_type objects_end,
			        MeshBase &mesh,
			        dofobject_accessor objects);

  /**
   * Distributes the global degrees of freedom, for dofs on
   * this processor.  In this format the local
   * degrees of freedom are in a contiguous block for each
   * variable in the system.
   * Starts at index next_free_dof, and increments it to
   * the post-final index.
   * If build_send_list is true, builds the send list.  If
   * false, clears and reserves the send list
   */  
  void distribute_local_dofs_var_major (unsigned int& next_free_dof,
				        MeshBase& mesh,
				        const bool build_send_list);
  
  /**
   * Distributes the global degrees of freedom, for dofs on
   * this processor.  In this format all the 
   * degrees of freedom at a node/element are in contiguous 
   * blocks.  Note in particular that the degrees of freedom
   * for a given variable are not in contiguous blocks, as 
   * in the case of \p distribute_local_dofs_var_major.
   * Starts at index next_free_dof, and increments it to
   * the post-final index.
   * If build_send_list is true, builds the send list.  If
   * false, clears and reserves the send list
   */  
  void distribute_local_dofs_node_major (unsigned int& next_free_dof,
				         MeshBase& mesh,
				         const bool build_send_list);
    
  /**
   * Adds entries to the \p _send_list vector corresponding to DoFs
   * on elements neighboring the current processor.
   */
  void add_neighbors_to_send_list(MeshBase& mesh);

  /**
   * Takes the \p _send_list vector (which may have duplicate entries)
   * and sorts it.  The duplicate entries are then removed, resulting in
   * a sorted \p _send_list with unique entries.
   */
  void sort_send_list ();

  
#if defined(ENABLE_AMR) || defined(ENABLE_PERIODIC)

  /**
   * Build the constraint matrix C associated with the element
   * degree of freedom indices elem_dofs. The optional parameter
   * \p called_recursively should be left at the default value
   * \p false.  This is used to handle the special case of
   * an element's degrees of freedom being constrained in terms
   * of other, local degrees of freedom.  The usual case is
   * for an elements DOFs to be constrained by some other,
   * external DOFs.   
   */
  void build_constraint_matrix (DenseMatrix<Number>& C,
				std::vector<unsigned int>& elem_dofs,
				const bool called_recursively=false) const;

  /**
   * Finds all the DOFS associated with the element DOFs elem_dofs.
   * This will account for off-element couplings via hanging nodes.
   */
  void find_connected_dofs (std::vector<unsigned int> &elem_dofs) const;
  
#endif


  /**
   * The finite element type for each variable.
   */
  std::vector<FEType*> _variable_types;

  /**
   * The number of the system we manage DOFs for.
   */
  const unsigned int _sys_number;
  
  /**
   * Additional matrices handled by this object.  These pointers do @e 
   * not handle the memory, instead, \p System, who
   * told \p DofMap about them, owns them.
   */
  std::vector<SparseMatrix<Number>* > _matrices;
  
  /**
   * First DOF index on processor \p p.
   */
  std::vector<unsigned int> _first_df;

  /**
   * Last DOF index (plus 1) on processor \p p.
   */
  std::vector<unsigned int> _end_df;

  /**
   * A list containing all the global DOF indicies that affect the
   * solution on my subdomain.
   */
  std::vector<unsigned int> _send_list;

  /**
   * The number of on-processor nonzeros in my portion of the
   * global matrix.
   */
  std::vector<unsigned int> _n_nz;
  
  /**
   * The number of off-processor nonzeros in my portion of the
   * global matrix.
   */
  std::vector<unsigned int> _n_oz;

  /**
   * Total number of degrees of freedom.
   */
  unsigned int _n_dfs;

#ifdef ENABLE_AMR

  /**
   * Total number of degrees of freedom on old dof objects
   */
  unsigned int _n_old_dfs;

#endif

#if defined(ENABLE_AMR) || defined(ENABLE_PERIODIC)

  /**
   * Data structure containing DOF constraints.  The ith
   * entry is the constraint matrix row for DOF i.
   */
  DofConstraints _dof_constraints;

#endif


#ifdef ENABLE_PERIODIC
  /**
   * Data structure containing periodic boundaries.  The ith
   * entry is the constraint matrix row for boundaryid i.
   */
  PeriodicBoundaries _periodic_boundaries;
#endif

  friend class SparsityPattern::Build;
};



// ------------------------------------------------------------
// Dof Map inline member functions
inline
DofMap::DofMap(const unsigned int number) :
  _dof_coupling(NULL),
  _sys_number(number),
//  _matrix(NULL),
  _n_dfs(0) 
#ifdef ENABLE_AMR
  , _n_old_dfs(0)
#endif
{
  _matrices.clear();
}


#if defined(ENABLE_AMR) || defined(ENABLE_PERIODIC)

inline
bool DofMap::is_constrained_dof (const unsigned int dof) const
{
  if (_dof_constraints.count(dof) != 0)
    return true;

  return false;
}

#endif


#ifdef ENABLE_PERIODIC

inline
bool DofMap::is_periodic_boundary (const unsigned int boundaryid) const
{
  if (_periodic_boundaries.count(boundaryid) != 0)
    return true;

  return false;
}

#endif


inline
unsigned int DofMap::sys_number() const
{
  return _sys_number;
}



// ------------------------------------------------------------
// SparsityPattern inline member functions
template<typename BidirectionalIterator>
inline
void SparsityPattern::sort_row (const BidirectionalIterator begin,
				BidirectionalIterator       middle,
				const BidirectionalIterator end)
{
  if ((begin == middle) || (middle == end)) return;

  libmesh_assert (std::distance (begin,  middle) > 0);
  libmesh_assert (std::distance (middle, end)    > 0);
  libmesh_assert (std::unique (begin,  middle) == middle);
  libmesh_assert (std::unique (middle, end)    == end);
  
  while (middle != end)
    {
      BidirectionalIterator
	b = middle,
	a = b-1;

      // Bubble-sort the middle value downward
      while (!(*a < *b)) // *a & *b are less-than comparable, so use <
	{
	  std::swap (*a, *b);

#if defined(__GNUC__) && (__GNUC__ < 4) && !defined(__INTEL_COMPILER)
	  /* Prohibit optimization at this point since gcc 3.3.5 seems
	     to have a bug.  */
	  SparsityPattern::_dummy_function();
#endif

	  if (a == begin) break;
	  
	  b=a;
	  --a;
	}
      
      ++middle;
    }

  // Assure the algorithm worked if we are in DEBUG mode
#ifdef DEBUG
  {
    // SGI STL extension!
    // libmesh_assert (std::is_sorted(begin,end));
    
    BidirectionalIterator
      prev  = begin,
      first = begin;

    for (++first; first != end; prev=first, ++first)
      if (*first < *prev)
	libmesh_assert(false);
  }  
#endif

  // Make sure the two ranges did not contain any common elements
  libmesh_assert (std::unique (begin, end) == end);
}


// ------------------------------------------------------------
// PeriodicBoundary inline member functions
#ifdef ENABLE_PERIODIC

inline
PeriodicBoundary *PeriodicBoundaries::boundary(unsigned int id)
{
  iterator i = this->find(id);
  if (i == this->end())
    return NULL;
  return &i->second;
}

#endif // ENABLE_PERIODIC


#endif // __dof_map_h__