#include #include #include #include "Debug.h" #include "DOFVector.h" #include "MacroElement.h" #include "VtkWriter.h" #include "ElementFileWriter.h" namespace AMDiS { namespace debug { #ifdef HAVE_PARALLEL_DOMAIN_AMDIS void writeLocalElementDofs(int rank, int elIdx, FiniteElemSpace *feSpace) { using boost::lexical_cast; if (MPI::COMM_WORLD.Get_rank() == rank) { DOFVector tmp(feSpace, "tmp"); colorDofVectorByLocalElementDofs(tmp, feSpace->getMesh(), elIdx); VtkWriter::writeFile(tmp, "tmp" + lexical_cast(elIdx) + ".vtu"); } } void writeDofMesh(int rank, DegreeOfFreedom dof, FiniteElemSpace *feSpace) { using boost::lexical_cast; if (MPI::COMM_WORLD.Get_rank() == rank) { DOFVector tmp(feSpace, "tmp"); tmp.set(0.0); tmp[dof] = 1.0; VtkWriter::writeFile(tmp, "dofmesh" + lexical_cast(rank) + ".vtu"); } } void writeMesh(FiniteElemSpace *feSpace, int rank, std::string filename) { using boost::lexical_cast; int myRank = MPI::COMM_WORLD.Get_rank(); if (rank == -1 || myRank == rank) { DOFVector tmp(feSpace, "tmp"); VtkWriter::writeFile(tmp, filename + lexical_cast(myRank) + ".vtu", false); } } #endif void colorDofVectorByLocalElementDofs(DOFVector& vec, Element *el) { // === Get local indices of the given element. === const BasisFunction *basisFcts = vec.getFESpace()->getBasisFcts(); int nBasisFcts = basisFcts->getNumber(); std::vector localDofs(nBasisFcts); basisFcts->getLocalIndices(el, vec.getFESpace()->getAdmin(), localDofs); // === Set the values of the dof vector. === vec.set(0.0); for (int i = 0; i < nBasisFcts; i++) vec[localDofs[i]] = static_cast(i); } bool colorDofVectorByLocalElementDofs(DOFVector& vec, Mesh *mesh, int elIndex) { FUNCNAME("colorDofVectorByLocalElementDofs()"); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { if (elInfo->getElement()->getIndex() == elIndex) { colorDofVectorByLocalElementDofs(vec, elInfo->getElement()); return true; } elInfo = stack.traverseNext(elInfo); } return false; } Element* getDofIndexElement(FiniteElemSpace *feSpace, DegreeOfFreedom dof) { const BasisFunction* basFcts = feSpace->getBasisFcts(); int nBasFcts = basFcts->getNumber(); std::vector dofVec(nBasFcts); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), dofVec); for (int i = 0; i < nBasFcts; i++) if (dofVec[i] == dof) return elInfo->getElement(); elInfo = stack.traverseNext(elInfo); } return NULL; } Element* getLevel0ParentElement(Mesh *mesh, Element *el) { TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { if (elInfo->getElement() == el) return elInfo->getMacroElement()->getElement(); elInfo = stack.traverseNext(elInfo); } return NULL; } void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof) { Element *el = getDofIndexElement(feSpace, dof); Element *parEl = getLevel0ParentElement(feSpace->getMesh(), el); std::cout << "DOF-INFO: dof = " << dof << " elidx = " << el->getIndex() << " pelidx = " << parEl->getIndex() << std::endl; TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { if (elInfo->getElement()->getIndex() == parEl->getIndex()) std::cout << "EL INFO TO " << parEl->getIndex() << ": " << elInfo->getType() << std::endl; elInfo = stack.traverseNext(elInfo); } } void printMatValuesStatistics(Matrix *mat) { std::map counter; int counter0 = 0; for (int i = 0; i < mat->getNumRows(); i++) { for (int j = 0; j < mat->getNumCols(); j++) { DOFMatrix *dofMat = (*mat)[i][j]; if (dofMat) { using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits = mtl::traits; typedef DOFMatrix::base_matrix_type base_matrix_type; traits::const_value::type value(dofMat->getBaseMatrix()); typedef traits::range_generator::type cursor_type; typedef traits::range_generator::type icursor_type; for (cursor_type cursor = begin(dofMat->getBaseMatrix()), cend = end(dofMat->getBaseMatrix()); cursor != cend; ++cursor) for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { if (value(*icursor) == 0.0) { counter0++; } else { int a = static_cast(floor(log10(fabs(value(*icursor))))); counter[a]++; } } } } } std::cout << "value = 0: " << counter0 << std::endl; for (std::map::iterator it = counter.begin(); it != counter.end(); ++it) std::cout << pow(10, it->first) << " <= values <= " << pow(10, it->first + 1) << ": " << it->second << std::endl; } void writeElementIndexMesh(FiniteElemSpace *feSpace, std::string filename) { std::map vec; TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(), -1, Mesh::CALL_LEAF_EL); while (elInfo) { int index = elInfo->getElement()->getIndex(); vec[index] = index; elInfo = stack.traverseNext(elInfo); } ElementFileWriter::writeFile(vec, feSpace, filename); } void writeMatlabMatrix(DOFMatrix &mat, std::string filename) { using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits = mtl::traits; typedef DOFMatrix::base_matrix_type Matrix; traits::row::type row(mat.getBaseMatrix()); traits::col::type col(mat.getBaseMatrix()); traits::const_value::type value(mat.getBaseMatrix()); typedef traits::range_generator::type cursor_type; typedef traits::range_generator::type icursor_type; std::ofstream out; out.open(filename.c_str()); for (cursor_type cursor = begin(mat.getBaseMatrix()), cend = end(mat.getBaseMatrix()); cursor != cend; ++cursor) for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) out << row(*icursor) + 1 << " " << col(*icursor) + 1 << " " << value(*icursor) << "\n"; out.close(); } void writeMatlabMatrix(Matrix &mat, std::string filename) { using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits = mtl::traits; typedef DOFMatrix::base_matrix_type Matrix; std::ofstream out; out.open(filename.c_str()); for (int i = 0; i < mat.getNumRows(); i++) { for (int j = 0; j < mat.getNumCols(); j++) { if (!mat[i][j]) continue; traits::row::type row(mat[i][j]->getBaseMatrix()); traits::col::type col(mat[i][j]->getBaseMatrix()); traits::const_value::type value(mat[i][j]->getBaseMatrix()); typedef traits::range_generator::type cursor_type; typedef traits::range_generator::type icursor_type; for (cursor_type cursor = begin(mat[i][j]->getBaseMatrix()), cend = end(mat[i][j]->getBaseMatrix()); cursor != cend; ++cursor) for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) out << i * num_rows(mat[i][j]->getBaseMatrix()) + row(*icursor) + 1 << " " << j * num_cols(mat[i][j]->getBaseMatrix()) + col(*icursor) + 1 << " " << value(*icursor) << "\n"; } } out.close(); } void writeMatlabVector(DOFVector &vec, std::string filename) { std::ofstream out; out.open(filename.c_str()); DOFIterator it(&vec, USED_DOFS); int counter = 1; for (it.reset(); !it.end(); ++it) out << counter++ << " " << *it << "\n"; out.close(); } void createSortedDofs(Mesh *mesh, ElementIdxToDofs &elMap) { FUNCNAME("Debug::dbgCreateElementMap()"); int dim = mesh->getDim(); elMap.clear(); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { Element *el = elInfo->getElement(); switch (dim) { case 2: sortDofs(el->getDOF(0), el->getDOF(1), el->getDOF(2), elMap[el->getIndex()]); break; case 3: sortDofs(el->getDOF(0), el->getDOF(1), el->getDOF(2), el->getDOF(3), elMap[el->getIndex()]); break; default: ERROR_EXIT("What is this?\n"); } elInfo = stack.traverseNext(elInfo); } } void testSortedDofs(Mesh *mesh, ElementIdxToDofs &elMap) { FUNCNAME("Debug::dbgTestElementMap()"); int dim = mesh->getDim(); int nVertex = Global::getGeo(VERTEX, dim); DofContainer vec(nVertex); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { Element *el = elInfo->getElement(); switch (dim) { case 2: sortDofs(el->getDOF(0), el->getDOF(1), el->getDOF(2), vec); break; case 3: sortDofs(el->getDOF(0), el->getDOF(1), el->getDOF(2), el->getDOF(3), vec); break; default: ERROR_EXIT("What is this?\n"); } for (int i = 0; i < nVertex; i++) { if (elMap[el->getIndex()][i] != vec[i]) { MSG("[DBG] Wrong new dof numeration in element = %d\n!", el->getIndex()); std::cout << "[DBG]: Old numeration was: "; for (int j = 0; j < nVertex; j++) std::cout << elMap[el->getIndex()][j] << " = " << *(elMap[el->getIndex()][j]) << " "; std::cout << std::endl; std::cout << "[DBG]: New numeration is: "; for (int j = 0; j < nVertex; j++) std::cout << vec[j] << " = " << *(vec[j]) << " "; std::cout << std::endl; ERROR_EXIT("WRONG NEW DOF NUMERATION!\n"); } } elInfo = stack.traverseNext(elInfo); } } void sortDofs(const DegreeOfFreedom* dof0, const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2, DofContainer &vec) { DofPtrSortFct dofPtrSort; vec.resize(3); vec[0] = dof0; vec[1] = dof1; vec[2] = dof2; sort(vec.begin(), vec.end(), dofPtrSort); } void sortDofs(const DegreeOfFreedom* dof0, const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2, const DegreeOfFreedom* dof3, DofContainer &vec) { DofPtrSortFct dofPtrSort; vec.resize(4); vec[0] = dof0; vec[1] = dof1; vec[2] = dof2; vec[3] = dof3; sort(vec.begin(), vec.end(), dofPtrSort); } } // namespace debug } // namespace AMDiS