-
Thomas Witkowski authoredThomas Witkowski authored
DOFMatrix.h 14.48 KiB
// ============================================================================
// == ==
// == AMDiS - Adaptive multidimensional simulations ==
// == ==
// == http://www.amdis-fem.org ==
// == ==
// ============================================================================
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.
/** \file DOFMatrix.h */
#ifndef AMDIS_DOFMATRIX_H
#define AMDIS_DOFMATRIX_H
#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 "DOFIndexed.h"
#include "Boundary.h"
#include "Serializable.h"
namespace AMDiS {
/** \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:
/// Type of scalars in the underlying matrix
typedef double value_type;
typedef unsigned size_type;
typedef mtl::matrix::parameters<mtl::row_major, mtl::index::c_index, mtl::non_fixed::dimensions, false, size_type> para;
/// Type of underlying matrix
typedef mtl::compressed2D<value_type, para> base_matrix_type;
/// Type of inserter for the base matrix;
typedef mtl::matrix::inserter<base_matrix_type, mtl::operations::update_plus<value_type> > inserter_type;
private:
/// 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();
/// Constructs a DOFMatrix with name n and the given row and olumn FeSpaces.
DOFMatrix(const FiniteElemSpace* rowFeSpace, const FiniteElemSpace* colFeSpace,
std::string n = "");
/// Copy-Constructor
DOFMatrix(const DOFMatrix& rhs);
/// Destructor
virtual ~DOFMatrix();
/// Assignment 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();
}
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;
}
/// 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(const ElementMatrix& elMat,
const BoundaryType *bound,
ElInfo* rowElInfo,
ElInfo* colElInfo);
/* \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);
/** \brief
* Finishes insertion. For compressed matrix types, this is where the
* compression happens.
*/
void finishInsertion()
{
FUNCNAME("DOFMatrix::finishInsertion()");
TEST_EXIT(inserter)("Inserter wasn't used or is already finished.\n");
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();
}
/// Returns \ref name
inline 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");
}
/// 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 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;
}
/// Calculates the average of non zero entries per row in matrix.
void calculateNnz()
{
nnzPerRow = 0;
if (num_rows(matrix) != 0)
nnzPerRow = int(double(matrix.nnz()) / num_rows(matrix) * 1.2);
if (nnzPerRow < 20)
nnzPerRow = 20;
}
/// Returns \ref nnzPerRow.
int getNnz()
{
return nnzPerRow;
}
/// Writes the matrix to an output stream.
void serialize(std::ostream &out);
/// Reads a matrix from an input stream.
void deserialize(std::istream &in);
///
int memsize();
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
void setRankDofs(std::map<DegreeOfFreedom, bool>& dofmap)
{
rankDofs = &dofmap;
}
#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;
/// Number of basis functions in the row fe space.
int nRow;
/// Number of basis functions in the col fe space.
int nCol;
/// Maps local row indices of an element to global matrix indices.
std::vector<DegreeOfFreedom> rowIndices;
/// Maps local col indices of an element to global matrix indices.
std::vector<DegreeOfFreedom> colIndices;
/* \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;
/* \brief
* Number of non zero entries per row (average). For instationary problems this
* information may be used in the next timestep to accelerate insertion of
* elemnts during assembling. */
int nnzPerRow;
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
/// Stores for the DOFs of the row FE spaces whether they are owned by the
/// rank or not. This is used to ensure that Dirichlet BC is handled
/// correctly in parallel computations.
std::map<DegreeOfFreedom, bool> *rankDofs;
#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