DOFMatrix.h

// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file DOFMatrix.h */

#ifndef AMDIS_DOFMATRIX_H
#define AMDIS_DOFMATRIX_H

// ============================================================================
// ===== includes =============================================================
// ============================================================================

#include <vector>
#include <set>
#include <memory>
#include <list>
#include <boost/numeric/mtl/mtl.hpp>

#include "AMDiS_fwd.h"
#include "Global.h"
#include "Flag.h"
#include "RCNeighbourList.h"
#include "DOFAdmin.h"
#include "DOFIterator.h"
#include "DOFIndexed.h"
#include "MemoryManager.h"
#include "Boundary.h"
#include "Serializable.h"

#ifdef HAVE_PARALLEL_AMDIS
#include "petscao.h"
#endif

namespace AMDiS {


  // ===========================================================================
  // ===== struct MatEntry =====================================================
  // ===========================================================================

  /** \brief
   * Represents one entry of a DOFMatrix
   */
  struct MatEntry
  {
    /** \brief
     * column index of entry (if >= 0; else unused)
     */
    DegreeOfFreedom col;

    /** \brief
     * matrix entry
     */
    double entry;
  };


  /** \brief
   * Is used to search for all entries of a matrix which column index is 
   * smaller than a given value.
   */
  class MatEntryValueLess : public std::unary_function<MatEntry, bool> {
  private:
    const double value_;
  public:
    MatEntryValueLess(const double& value) 
      : value_(value) 
    {}

    bool operator()(const MatEntry& m) {
      return (fabs(m.entry) < value_);
    }      
  };


  /** \brief
   * Is used to search for all entries of a matrix which value is smaller 
   * than a given value.
   */
  class MatEntryColLess : public std::unary_function<MatEntry, bool> {
  private:
    const int col_;
  public:
    MatEntryColLess(const int& col) 
      : col_(col) 
    {}

    bool operator()(const MatEntry& m) {
      return (fabs(m.col) < col_);
    }      
  };


  /** \brief
   * This function is required if two MatEntries are compared by their col
   * entry (for example when sorting a vector of MatEntries).
   */
  struct CmpMatEntryCol : public std::binary_function<MatEntry, MatEntry, bool> {
    bool operator()(const MatEntry& m1, const MatEntry m2) const {
      return m1.col < m2.col;
    }
  };
  

  /** \brief
   * This function compares two matrix entries by their values. It returns true,
   * if the value of m2 is greater than the value of m1.
   */
  struct CmpMatEntryValue : public std::binary_function<MatEntry, MatEntry, bool> {
    bool operator()(const MatEntry& m1, const MatEntry m2) const {
      return m1.entry < m2.entry;
    }
  };


  /** \brief
   * This function compares two matrix entries by their values. It returns true,
   * if the value of m2 is greater than the value of m1.
   */
  struct CmpMatEntryAbsValueLess : public std::binary_function<MatEntry, MatEntry, bool> {
    bool operator()(const MatEntry& m1, const MatEntry m2) const {
      return fabs(m1.entry) < fabs(m2.entry);
    }
  };
  
  /** \brief
   * This function compares two matrix entries by their values. It returns true,
   * if the value of m1 is greater than the value of m2.
   */
  struct CmpMatEntryAbsValueGreater : public std::binary_function<MatEntry, MatEntry, bool> {
    bool operator()(const MatEntry& m1, const MatEntry m2) const {
      return fabs(m1.entry) > fabs(m2.entry);
    }
  };

  // ============================================================================
  // ===== class DOFMatrix ======================================================
  // ============================================================================

  /** \ingroup DOFAdministration
   * \brief
   * A DOFMatrix is a sparse matrix representation for matrices that work
   * on DOFVectors. The underlying matrix type 
   * is a CRS matrix of double.
   */
  class DOFMatrix : public DOFIndexed<bool>,
                    public Serializable
  {
  public:
    MEMORY_MANAGED(DOFMatrix);

    /// Type of scalars in the underlying matrix
    typedef double                          value_type;

    /// Type of underlying matrix
    typedef mtl::compressed2D<value_type>       base_matrix_type;

    /// Type of inserter for the base matrix;
    typedef mtl::matrix::inserter<base_matrix_type>  inserter_type;

  private:
    /** \ingroup DOFAdministration
     * \brief
     * Alias for DOFIterator< std::vector<MatEntry<T> > >. So one can access an
     * iterator working on a DOFMatrix via DOFMatrix::Iterator
     */
    class Iterator : public DOFIterator< std::vector<MatEntry> > {
    public:
      Iterator(DOFIndexed< std::vector<MatEntry> > *c, 
	       DOFIteratorType type)
	: DOFIterator< std::vector<MatEntry> >(c, type)
      {}
    };
    
    /** \brief
     * A MatrixRow represents one row of the sparse matrix. It is realized by
     * a STL vector of MatEntry<T> objects
     */
    typedef std::vector<MatEntry> MatrixRow;
    
    /** \brief
     * Symbolic constant for an unused matrix entry
     */
    static const int UNUSED_ENTRY = -1;
    
    /** \brief
     * Symbolic constant for an unused entry which is not followed by any
     * used entry in this row
     */
    static const int NO_MORE_ENTRIES = -2;

  public:    
    DOFMatrix();

    /** \brief
     * Constructs a DOFMatrix with name n and the given row and 
     * column FESpaces.
     */
    DOFMatrix(const FiniteElemSpace* rowFESpace,
	      const FiniteElemSpace* colFESpace,
	      std::string n = "");

    /** \brief
     * Copy-Constructor
     */
    DOFMatrix(const DOFMatrix& rhs);


    /** \brief
     * Destructor
     */
    virtual ~DOFMatrix();
  
    /** \brief
     * operator=
     */
    DOFMatrix& operator=(const DOFMatrix& rhs);

    void copy(const DOFMatrix& rhs);

    /// Access underlying matrix directly
    base_matrix_type& getBaseMatrix()
    {
	return matrix;
    }

    /// Access underlying matrix directly (const)
    const base_matrix_type& getBaseMatrix() const
    {
	return matrix;
    }

    // Only to get rid of the abstract functions, I hope they are not used
    std::vector<bool>::iterator begin() {ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v; return v.begin();}
    std::vector<bool>::iterator end() {ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v; return v.end();}
#if 1
    
    bool dummy; // Must be deleted later
    bool& operator[](int i) {ERROR_EXIT("Shouldn't be used, only fake."); return dummy;}
    const bool& operator[](int i) const {ERROR_EXIT("Shouldn't be used, only fake."); return dummy;}
#endif
 
    /// DOFMatrix does not need to be compressed before assembling, when using MTL4.
    void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF) {}

    /// Implements DOFIndexedBase::freeDOFContent()
    virtual void freeDOFContent(int index);

    /// Returns \ref coupleMatrix.
    inline bool isCoupleMatrix() { 
      return coupleMatrix; 
    }

    /// Returns \ref coupleMatrix.
    inline void setCoupleMatrix(bool c) { 
      coupleMatrix = c; 
    }

    /// a * x + y
    void axpy(double a, const DOFMatrix& x, const DOFMatrix& y);

    /// Multiplication with a scalar.
    void scal(double s);

    /** \brief
     * Adds an operator to the DOFMatrix. A factor, that is multipled
     * to the operator, and a multilier factor for the estimator may be
     * also given.
     */
    void addOperator(Operator *op, double* factor = NULL, double* estFactor = NULL);

    inline std::vector<double*>::iterator getOperatorFactorBegin() {
      return operatorFactor.begin();
    }

    inline std::vector<double*>::iterator getOperatorFactorEnd() {
      return operatorFactor.end();
    }

    inline std::vector<double*>::iterator getOperatorEstFactorBegin() {
      return operatorEstFactor.begin();
    }

    inline std::vector<double*>::iterator getOperatorEstFactorEnd() {
      return operatorEstFactor.end();
    }

    inline std::vector<Operator*>::iterator getOperatorsBegin() {
      return operators.begin();
    }

    inline std::vector<Operator*>::iterator getOperatorsEnd() {
      return operators.end();
    }

    Flag getAssembleFlag();

    /** \brief
     * Updates the matrix matrix by traversing the underlying mesh and assembling
     * the element contributions into the matrix. Information about the 
     * computation of element matrices and connection of local and global DOFs is
     * stored in minfo; the flags for the mesh traversal are stored at 
     * minfo->fill_flags which specifies the elements to be visited and 
     * information that should be present on the elements for the calculation of 
     * the element matrices and boundary information (if minfo->boundBas is not
     * NULL). On the elements, information about the row DOFs is accessed by 
     * minfo->rowBas using info->row_admin; this vector is also used for the 
     * column DOFs if minfo->nCol is less or equal zero, or minfo->col_admin or 
     * minfo->colBas is a NULL pointer; if row and column DOFs are the same, the 
     * boundary type of the DOFs is accessed by minfo->boundBas if 
     * minfo->boundBas is not a NULL pointer; then the element matrix is 
     * computed by minfo->fillElementMatrix(el info, minfo->myFill); these 
     * contributions, multiplied by minfo->factor, are eventually added to matrix
     * by a call of addElementMatrix() with all information about row and column 
     * DOFs, the element matrix, and boundary types, if available;
     * updateMatrix() only adds element contributions; this makes several calls 
     * for the assemblage of one matrix possible; before the first call, the 
     * matrix should be cleared by calling clear dof matrix().
     */
  
    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound);

    void assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
		  Operator *op);

    void assemble(double factor, 
		  ElInfo *rowElInfo, ElInfo *colElInfo,
		  ElInfo *smallElInfo, ElInfo *largeElInfo,
		  const BoundaryType *bound,
		  Operator *op = NULL);

    void assemble2(double factor, 
		   ElInfo *mainElInfo, ElInfo *auxElInfo,
		   ElInfo *smallElInfo, ElInfo *largeElInfo,
		   const BoundaryType *bound, 
		   Operator *op = NULL);

    /// Adds an element matrix to \ref matrix
    void addElementMatrix(double sign, 
			  const ElementMatrix& elMat, 
			  const BoundaryType *bound,
			  bool add = true);

    /* \brief
     * That function must be called after the matrix assembling has been finished. 
     * This makes it possible to start some cleanup or matrix data compressing 
     * procedures.
     */
    void finishAssembling();

    /** \brief
     * Enable insertion for assembly. You can optionally give an upper limit for
     * entries per row (per column for CCS matrices).  Choosing this parameter
     * too small can induce perceivable overhead for compressed matrices.  Thus,
     * it's better to choose a bit too large than too small.
     */
    void startInsertion(int nnz_per_row= 10)
    {
      if (inserter) 
	delete inserter, inserter= 0; 

      inserter= new inserter_type(matrix, nnz_per_row);
    }

    /** \brief
     *  Finishes insertion. For compressed matrix types, this is where the 
     *  compression happens.
     */
    void finishInsertion()
    {
      TEST_EXIT(inserter)("Inserter wasn't used or is already finished.");
      
      delete inserter;
      inserter= 0;
    }

    /** \brief
     * Returns whether restriction should be performed after coarsening
     * (false by default)
     */
    virtual bool coarseRestrict() {
      return false;
    }

    /// Returns const \ref rowFESpace
    const FiniteElemSpace* getRowFESpace() const { 
      return rowFESpace; 
    }

    /// Returns const \ref colFESpace
    const FiniteElemSpace* getColFESpace() const { 
      return colFESpace; 
    }

    /// Returns const \ref rowFESpace
    const FiniteElemSpace* getFESpace() const { 
      return rowFESpace; 
    }

    /// Returns number of rows (\ref matrix.size())
    inline int getSize() const { 
      return num_rows(matrix);
    }

    /** \brief
     * Returns the number of used rows (equal to number of used DOFs in
     * the row FE space).
     */
    inline int getUsedSize() const {
      return rowFESpace->getAdmin()->getUsedSize();
    }

    // Only fake, shouldn't be called
    /** \brief
     * Returns number of cols. For that, the function iteratos over all
     * rows and searchs for the entry with the highest col number.
     */
    int getNumCols() const;

    /// Returns \ref name
    inline const std::string& getName() const { 
      return name; 
    }

    /// Resizes \ref matrix to n rows
    inline void resize(int n) { 
      TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n");
      //      matrix.change_dim(n, n);
      //WARNING("Obsolete function without effect -- Peter\n");
    }

    /// Returns value at logical indices a,b
    double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;

    /// Changes col at logical indices a,b to c 
    void changeColOfEntry(DegreeOfFreedom a, DegreeOfFreedom b, DegreeOfFreedom c);

    /** \brief
     * Creates an entry with logical indices irow, icol if there is no entry
     * yet. Than sign * entry is added to the value at this logical indices
     */
    void addSparseDOFEntry(double sign,
			   int irow, int jcol, double entry,
			   bool add = true);

    void addMatEntry(int row, MatEntry entry);

    void addMatEntry(int row, int DegreeOfFreedom, double value);

    void addRow(std::vector<MatEntry> row);

    void removeRowsWithDBC(std::set<int> *rows);

    /// Prints \ref matrix to stdout
    void print() const;

    /// Removes all matrix entries
    void clear()
    {
	set_to_zero(matrix);
    }

    /// Test whether \ref matrix is symmetric. Exits if not.
    void test();

    bool symmetric();

    inline std::vector<Operator*> getOperators() { 
      return operators; 
    }
    
    inline std::vector<double*> getOperatorFactor() { 
      return operatorFactor; 
    }

    inline std::vector<double*> getOperatorEstFactor() { 
      return operatorEstFactor; 
    }

    inline BoundaryManager* getBoundaryManager() const { 
      return boundaryManager; 
    }

    /// Returns a pointer to \ref applyDBCs.
    std::set<int>* getApplyDBCs() {
      return &applyDBCs;
    }

    inline void setBoundaryManager(BoundaryManager *bm) {
      boundaryManager = bm;
    }

    void createPictureFile(const char* filename, int dim);
   
  private:
    template <typename T>
    void s_write(::std::ostream &out, const T& value)
    {
      out.write(reinterpret_cast<const char*>(&value), sizeof value);
    }

  public:    
    void serialize(::std::ostream &out) 
    {
      using namespace mtl; 

      typedef traits::range_generator<tag::major, base_matrix_type>::type c_type;
      typedef traits::range_generator<tag::nz, c_type>::type              ic_type;

      typedef Collection<base_matrix_type>::size_type                     size_type;
      typedef Collection<base_matrix_type>::value_type                    value_type;
      
      traits::row<base_matrix_type>::type                                 row(matrix); 
      traits::col<base_matrix_type>::type                                 col(matrix);
      traits::const_value<base_matrix_type>::type                         value(matrix); 

      size_type rows= num_rows(matrix), cols= num_cols(matrix), total= matrix.nnz();
      s_write(out, rows);
      s_write(out, cols);
      s_write(out, total);

      for (c_type cursor(mtl::begin<tag::major>(matrix)), cend(mtl::end<tag::major>(matrix)); 
	   cursor != cend; ++cursor)
	for (ic_type icursor(mtl::begin<tag::nz>(cursor)), icend(mtl::end<tag::nz>(cursor)); 
	     icursor != icend; ++icursor) {
	  size_type   my_row= row(*icursor), my_col= col(*icursor);
	  value_type  my_value= value(*icursor);
	  s_write(out, my_row);
	  s_write(out, my_col);
	  s_write(out, my_value);
	}
    }

  private:
    template <typename T>
    void s_read(::std::istream &in, T& value)
    {
      in.read(reinterpret_cast<char*>(&value), sizeof value);
    }

  public:    
    void deserialize(::std::istream &in) 
    {
      using namespace mtl;

      typedef Collection<base_matrix_type>::size_type                     size_type;
      typedef Collection<base_matrix_type>::value_type                    value_type;
      
      size_type rows, cols, total;
      s_read(in, rows);
      s_read(in, cols);
      s_read(in, total);

      // Prepare matrix insertion
      clear();
      // matrix.change_dim(rows, cols) // do we want this?
      inserter_type ins(matrix);

      for (size_type i= 0; i < total; ++i) {
	size_type   my_row, my_col;
	value_type  my_value;
	s_read(in, my_row);
	s_read(in, my_col);
	s_read(in, my_value);
	ins(my_row, my_col) << my_value;
      }
    }

    ///
    int memsize();

#ifdef HAVE_PARALLEL_AMDIS
    /// Sets the petsc application ordering object to map dof indices.
    void useApplicationOrdering(AO *ao) {
      applicationOrdering = ao;
    }
#endif

  protected:
    /** \brief
     * Pointer to a FiniteElemSpace with information about corresponding row DOFs
     * and basis functions
     */
    const FiniteElemSpace *rowFESpace;

    /** \brief
     * Pointer to a FiniteElemSpace with information about corresponding 
     * column DOFs and basis functions
     */
    const FiniteElemSpace *colFESpace;

    /// Name of the DOFMatrix
    std::string name;

    /// Sparse matrix, type is a template parameter by default compressed2D<double>
    base_matrix_type matrix;

    /// Used while mesh traversal
    static DOFMatrix *traversePtr;
  
    /** \brief
     * Pointers to all operators of the equation systems. Are used in the
     * assembling process.
     */
    std::vector<Operator*> operators;
    
    /** \brief
     * Defines for each operator a factor which is used to scal the element
     * matrix after the assembling process of the operator.
     */
    std::vector<double*> operatorFactor;

    ///
    std::vector<double*> operatorEstFactor;

    ///
    BoundaryManager *boundaryManager;

    /** \brief
     * If false, the matrix is a diagonal matrix within a matrix of DOF matrices.
     * Otherwise the value is true, and the matrix is an off-diagonal matrix.
     */
    bool coupleMatrix;

    /// Temporary variable used in assemble()
    ElementMatrix *elementMatrix;

    /* \brief
     * A set of row indices. When assembling the DOF matrix, all rows, that
     * correspond to a dof at a dirichlet boundary, are ignored and the row is
     * left blank. After assembling, the diagonal element of the matrix must be
     * set to 1. The indices of all rows, where this should be done, are stored
     * in this set.
     */    
    std::set<int> applyDBCs;

#ifdef HAVE_PARALLEL_AMDIS
    /// Petsc application ordering to map dof indices.
    AO *applicationOrdering;
#endif

    /// Inserter object: implemented as pointer, allocated and deallocated as needed
    inserter_type *inserter;
      
    friend class DOFAdmin;
    friend class DOFVector<double>;
    friend class DOFVector<unsigned char>;
    friend class DOFVector<int>;
    friend class DOFVector<WorldVector<double> >;
  };

}

#endif  // AMDIS_DOFMATRIX_H