#include <algorithm> #include <png.h> #include <boost/numeric/mtl/mtl.hpp> #include "DOFMatrix.h" #include "QPsiPhi.h" #include "BasisFunction.h" #include "Boundary.h" #include "DOFAdmin.h" #include "ElInfo.h" #include "FiniteElemSpace.h" #include "Mesh.h" #include "DOFVector.h" #include "Operator.h" #include "BoundaryCondition.h" #include "BoundaryManager.h" #include "Assembler.h" #include "Utilities.h" namespace AMDiS { using namespace mtl; DOFMatrix *DOFMatrix::traversePtr = NULL; DOFMatrix::DOFMatrix() : rowFESpace(NULL), colFESpace(NULL), elementMatrix(3, 3), nRow(0), nCol(0), inserter(NULL) {} DOFMatrix::DOFMatrix(const FiniteElemSpace* rowFESpace_, const FiniteElemSpace* colFESpace_, std::string name_) : rowFESpace(rowFESpace_), colFESpace(colFESpace_), name(name_), coupleMatrix(false), inserter(NULL) { TEST_EXIT(rowFESpace)("no rowFESpace\n"); if (!colFESpace) colFESpace = rowFESpace; if (rowFESpace && rowFESpace->getAdmin()) (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->addDOFIndexed(this); boundaryManager = new BoundaryManager(rowFESpace_); nRow = rowFESpace->getBasisFcts()->getNumber(); nCol = colFESpace->getBasisFcts()->getNumber(); elementMatrix.change_dim(nRow, nCol); rowIndices.resize(nRow); colIndices.resize(nCol); applyDBCs.clear(); #ifdef HAVE_PARALLEL_DOMAIN_AMDIS applicationOrdering = NULL; #endif } DOFMatrix::DOFMatrix(const DOFMatrix& rhs) : name(rhs.name + "copy") { *this = rhs; if (rowFESpace && rowFESpace->getAdmin()) (const_cast<DOFAdmin*>( rowFESpace->getAdmin()))->addDOFIndexed(this); TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion"); inserter= 0; } DOFMatrix::~DOFMatrix() { FUNCNAME("DOFMatrix::~DOFMatrix()"); if (rowFESpace && rowFESpace->getAdmin()) (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this); if (boundaryManager) delete boundaryManager; if (inserter) delete inserter; } void DOFMatrix::print() const { FUNCNAME("DOFMatrix::print()"); using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits= mtl::traits; typedef base_matrix_type Matrix; traits::row<Matrix>::type row(matrix); traits::col<Matrix>::type col(matrix); traits::const_value<Matrix>::type value(matrix); typedef traits::range_generator<major, Matrix>::type cursor_type; typedef traits::range_generator<nz, cursor_type>::type icursor_type; std::cout.precision(10); for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor) { for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) if (value(*icursor) != 0.0) std::cout << "(" << row(*icursor) << "," << col(*icursor) << "," << value(*icursor) << ") "; std::cout << "\n"; } } bool DOFMatrix::symmetric() { FUNCNAME("DOFMatrix::symmetric()"); double tol = 1e-5; using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits= mtl::traits; typedef base_matrix_type Matrix; traits::row<Matrix>::type row(matrix); traits::col<Matrix>::type col(matrix); traits::const_value<Matrix>::type value(matrix); typedef traits::range_generator<major, Matrix>::type cursor_type; typedef traits::range_generator<nz, cursor_type>::type icursor_type; for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor) for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) // Compare each non-zero entry with its transposed if (abs(value(*icursor) - matrix[col(*icursor)][row(*icursor)]) > tol) return false; return true; } void DOFMatrix::test() { FUNCNAME("DOFMatrix::test()"); int non_symmetric = !symmetric(); if (non_symmetric) { MSG("matrix `%s' not symmetric.\n", name.data()); } else { MSG("matrix `%s' is symmetric.\n", name.data()); } } DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs) { rowFESpace = rhs.rowFESpace; colFESpace = rhs.colFESpace; operators = rhs.operators; operatorFactor = rhs.operatorFactor; coupleMatrix = rhs.coupleMatrix; /// The matrix values may only be copyed, if there is no active insertion. if (rhs.inserter == 0 && inserter == 0) matrix = rhs.matrix; if (rhs.boundaryManager) { boundaryManager = new BoundaryManager(*rhs.boundaryManager); } else { boundaryManager = NULL; } nRow = rhs.nRow; nCol = rhs.nCol; elementMatrix.change_dim(nRow, nCol); return *this; } void DOFMatrix::addElementMatrix(double sign, const ElementMatrix& elMat, const BoundaryType *bound, ElInfo* rowElInfo, ElInfo* colElInfo, bool add) { FUNCNAME("DOFMatrix::addElementMatrix()"); TEST_EXIT_DBG(inserter)("DOFMatrix is not in insertion mode"); inserter_type &ins= *inserter; // === Get indices mapping from local to global matrix indices. === rowFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(), rowFESpace->getAdmin(), &rowIndices); if (rowFESpace == colFESpace) { colIndices = rowIndices; } else { if (colElInfo) { colFESpace->getBasisFcts()->getLocalIndicesVec(colElInfo->getElement(), colFESpace->getAdmin(), &colIndices); } else { // If there is no colElInfo pointer, use rowElInfo the get the indices. colFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(), colFESpace->getAdmin(), &colIndices); } } // === Add element matrix to the global matrix using the indices mapping. === for (int i = 0; i < nRow; i++) { // for all rows of element matrix DegreeOfFreedom row = rowIndices[i]; BoundaryCondition *condition = bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL; if (condition && condition->isDirichlet()) { if (condition->applyBoundaryCondition()) applyDBCs.insert(static_cast<int>(row)); } else for (int j = 0; j < nCol; j++) { // for all columns DegreeOfFreedom col = colIndices[j]; double entry = elMat[i][j]; if (add) ins[row][col]+= sign * entry; else ins[row][col]= sign * entry; } } } double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const { return matrix[a][b]; } void DOFMatrix::freeDOFContent(int index) {} void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound) { FUNCNAME("DOFMatrix::assemble()"); set_to_zero(elementMatrix); std::vector<Operator*>::iterator it = operators.begin(); std::vector<double*>::iterator factorIt = operatorFactor.begin(); for (; it != operators.end(); ++it, ++factorIt) if ((*it)->getNeedDualTraverse() == false) (*it)->getElementMatrix(elInfo, elementMatrix, *factorIt ? **factorIt : 1.0); addElementMatrix(factor, elementMatrix, bound, elInfo, NULL); } void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound, Operator *op) { FUNCNAME("DOFMatrix::assemble()"); TEST_EXIT_DBG(op)("No operator!\n"); set_to_zero(elementMatrix); op->getElementMatrix(elInfo, elementMatrix, factor); addElementMatrix(factor, elementMatrix, bound, elInfo, NULL); } void DOFMatrix::assemble(double factor, ElInfo *rowElInfo, ElInfo *colElInfo, ElInfo *smallElInfo, ElInfo *largeElInfo, const BoundaryType *bound, Operator *op) { FUNCNAME("DOFMatrix::assemble()"); if (!op && operators.size() == 0) return; set_to_zero(elementMatrix); if (op) { op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, elementMatrix); } else { std::vector<Operator*>::iterator it = operators.begin(); std::vector<double*>::iterator factorIt = operatorFactor.begin(); for (; it != operators.end(); ++it, ++factorIt) { (*it)->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo, elementMatrix, *factorIt ? **factorIt : 1.0); } } addElementMatrix(factor, elementMatrix, bound, rowElInfo, colElInfo); } void DOFMatrix::assemble2(double factor, ElInfo *mainElInfo, ElInfo *auxElInfo, ElInfo *smallElInfo, ElInfo *largeElInfo, const BoundaryType *bound, Operator *op) { FUNCNAME("DOFMatrix::assemble2()"); if (!op && operators.size() == 0) { return; } set_to_zero(elementMatrix); if (op) { ERROR_EXIT("TODO"); // op->getElementMatrix(rowElInfo, colElInfo, // smallElInfo, largeElInfo, // elementMatrix); } else { std::vector<Operator*>::iterator it; std::vector<double*>::iterator factorIt; for(it = operators.begin(), factorIt = operatorFactor.begin(); it != operators.end(); ++it, ++factorIt) { if ((*it)->getNeedDualTraverse()) { (*it)->getElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, elementMatrix, *factorIt ? **factorIt : 1.0); } } } addElementMatrix(factor, elementMatrix, bound, mainElInfo, NULL); } void DOFMatrix::finishAssembling() { // call the operatos cleanup procedures for (std::vector<Operator*>::iterator it = operators.begin(); it != operators.end(); ++it) (*it)->finishAssembling(); } // Should work as before Flag DOFMatrix::getAssembleFlag() { Flag fillFlag(0); for (std::vector<Operator*>::iterator op = operators.begin(); op != operators.end(); ++op) fillFlag |= (*op)->getFillFlag(); return fillFlag; } void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y) { matrix+= a * x.matrix + y.matrix; } void DOFMatrix::scal(double b) { matrix*= b; } void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor) { operators.push_back(op); operatorFactor.push_back(factor); operatorEstFactor.push_back(estFactor); } void DOFMatrix::copy(const DOFMatrix& rhs) { matrix= rhs.matrix; } void DOFMatrix::removeRowsWithDBC(std::set<int> *rows) { inserter_type &ins= *inserter; for (std::set<int>::iterator it = rows->begin(); it != rows->end(); ++it) ins[*it][*it] = 1.0; rows->clear(); } int DOFMatrix::memsize() { return (num_rows(matrix) + matrix.nnz()) * sizeof(base_matrix_type::size_type) + matrix.nnz() * sizeof(base_matrix_type::value_type); } }