numeric_vector.h

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

  /**
   * \f$ U=v \f$ where v is a DenseVector<T> 
   * and you want to specify WHERE to insert it
   */
  virtual void insert (const std::vector<T>& v,
		       const std::vector<unsigned int>& dof_indices) = 0;

  /**
   * \f$U=V\f$, where U and V are type 
   * NumericVector<T> and you
   * want to specify WHERE to insert
   * the NumericVector<T> V 
   */
  virtual void insert (const NumericVector<T>& V,
		       const std::vector<unsigned int>& dof_indices) = 0;
      
  /**
   * \f$ U+=V \f$ where U and V are type 
   * DenseVector<T> and you
   * want to specify WHERE to insert
   * the DenseVector<T> V 
   */
  virtual void insert (const DenseVector<T>& V,
		       const std::vector<unsigned int>& dof_indices) = 0;
    
  /**
   * Scale each element of the
   * vector by the given factor.
   */
  virtual void scale (const T factor) = 0;

  /**
   * Computes the dot product, p = U.V
   */
  virtual T dot(const NumericVector<T>&) const = 0;
  
  /**
   * Creates a copy of the global vector in the
   * local vector \p v_local.
   */
  virtual void localize (std::vector<T>& v_local) const = 0;

  /**
   * Same, but fills a \p NumericVector<T> instead of
   * a \p std::vector.
   */
  virtual void localize (NumericVector<T>& v_local) const = 0;

  /**
   * Creates a local vector \p v_local containing
   * only information relevant to this processor, as
   * defined by the \p send_list.
   */
  virtual void localize (NumericVector<T>& v_local,
			 const std::vector<unsigned int>& send_list) const = 0;
  
  /**
   * Updates a local vector with selected values from neighboring
   * processors, as defined by \p send_list.
   */
  virtual void localize (const unsigned int first_local_idx,
			 const unsigned int last_local_idx,
			 const std::vector<unsigned int>& send_list) = 0;

  /**
   * Creates a local copy of the global vector in
   * \p v_local only on processor \p proc_id.  By
   * default the data is sent to processor 0.  This method
   * is useful for outputting data from one processor.
   */
  virtual void localize_to_one (std::vector<T>& v_local,
				const unsigned int proc_id=0) const = 0;
    
  /**
   * @returns \p -1 when \p this is equivalent to \p other_vector,
   * up to the given \p threshold.  When differences occur,
   * the return value contains the first index where
   * the difference exceeded the threshold.  When
   * no threshold is given, the \p libMesh \p TOLERANCE
   * is used.
   */
  virtual int compare (const NumericVector<T> &other_vector,
		       const Real threshold = TOLERANCE) const;

  /**
   * Prints the local contents of the vector to the screen.
   */
  virtual void print(std::ostream& os=std::cout) const;

  /**
   * Prints the global contents of the vector to the screen.
   */
  virtual void print_global(std::ostream& os=std::cout) const;

  /**
   * Same as above but allows you to use stream syntax.
   */
  friend std::ostream& operator << (std::ostream& os, const NumericVector<T>& v)
  {
    v.print_global(os);
    return os;
  }
  
  /**
   * Print the contents of the matrix in Matlab's
   * sparse matrix format. Optionally prints the
   * matrix to the file named \p name.  If \p name
   * is not specified it is dumped to the screen.
   */
  virtual void print_matlab(const std::string name="NULL") const
  {
    std::cerr << "ERROR: Not Implemented in base class yet!" << std::endl;
    std::cerr << "ERROR writing MATLAB file " << name << std::endl;
    libmesh_error();
  }

  /**
   * Creates the subvector "subvector" from the indices in the
   * "rows" array.  Similar to the create_submatrix routine for
   * the SparseMatrix class, it is currently only implemented for
   * PetscVectors.
   */
  virtual void create_subvector(NumericVector<T>& ,
				const std::vector<unsigned int>& ) const
  {
    std::cerr << "ERROR: Not Implemented in base class yet!" << std::endl;
    libmesh_error();
  }
    
protected:
  
  /**
   * Flag to see if the Numeric
   * assemble routines have been called yet
   */
  bool _is_closed;
  
  /**
   * Flag to tell if init 
   * has been called yet
   */
  bool _is_initialized;
};


/*----------------------- Inline functions ----------------------------------*/



template <typename T>
inline
NumericVector<T>::NumericVector () :
  _is_closed(false),
  _is_initialized(false)
{}



template <typename T>
inline
NumericVector<T>::NumericVector (const unsigned int n) :
  _is_closed(false),
  _is_initialized(false)
{
  init(n, n, false);
}



template <typename T>
inline
NumericVector<T>::NumericVector (const unsigned int n,
				 const unsigned int n_local) :
  _is_closed(false),
  _is_initialized(false)
{
  init(n, n_local, false);
}



template <typename T>
inline
NumericVector<T>::~NumericVector ()
{
  clear ();
}



// These should be pure virtual, not bugs waiting to happen - RHS
/*
template <typename T>
inline
NumericVector<T> & NumericVector<T>::operator= (const T) 
{
  //  libmesh_error();

  return *this;
}



template <typename T>
inline
NumericVector<T> & NumericVector<T>::operator= (const NumericVector<T>&) 
{
  //  libmesh_error();

  return *this;
}



template <typename T>
inline
NumericVector<T> & NumericVector<T>::operator= (const std::vector<T>&) 
{
  //  libmesh_error();

  return *this;
}
*/



template <typename T>
inline
void NumericVector<T>::clear ()
{
  _is_closed      = false;
  _is_initialized = false;
}



// Full specialization of the print() member for complex
// variables.  This must precede the non-specialized
// version, at least according to icc v7.1
template <>
inline
void NumericVector<Complex>::print(std::ostream& os) const
{
  libmesh_assert (this->initialized());
  os << "Size\tglobal =  " << this->size()
     << "\t\tlocal =  " << this->local_size() << std::endl;
  
  // std::complex<>::operator<<() is defined, but use this form
  os << "#\tReal part\t\tImaginary part" << std::endl;
  for (unsigned int i=this->first_local_index(); i<this->last_local_index(); i++)
    os << i << "\t" 
       << (*this)(i).real() << "\t\t" 
       << (*this)(i).imag() << std::endl;
}



template <typename T>
inline
void NumericVector<T>::print(std::ostream& os) const
{
  libmesh_assert (this->initialized());
  os << "Size\tglobal =  " << this->size()
     << "\t\tlocal =  " << this->local_size() << std::endl;

  os << "#\tValue" << std::endl;
  for (unsigned int i=this->first_local_index(); i<this->last_local_index(); i++)
    os << i << "\t" << (*this)(i) << std::endl;
}



template <>
inline
void NumericVector<Complex>::print_global(std::ostream& os) const
{
  libmesh_assert (this->initialized());

  std::vector<Complex> v(this->size());
  this->localize(v);

  // Right now we only want one copy of the output
  if (libMesh::processor_id())
    return;
  
  os << "Size\tglobal =  " << this->size() << std::endl;
  os << "#\tReal part\t\tImaginary part" << std::endl;
  for (unsigned int i=0; i!=v.size(); i++)
    os << i << "\t" 
       << v[i].real() << "\t\t" 
       << v[i].imag() << std::endl;
}


template <typename T>
inline
void NumericVector<T>::print_global(std::ostream& os) const
{
  libmesh_assert (this->initialized());

  std::vector<T> v(this->size());
  this->localize(v);

  // Right now we only want one copy of the output
  if (libMesh::processor_id())
    return;

  os << "Size\tglobal =  " << this->size() << std::endl;
  os << "#\tValue" << std::endl;
  for (unsigned int i=0; i!=v.size(); i++)
    os << i << "\t" << v[i] << std::endl;
}



#endif  // #ifdef __numeric_vector_h__