diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc index ed2338e5d3893f120cc444603c8feb82ec296f8b..27b3a03ac391903b39e781c82f7938231943f117 100644 --- a/AMDiS/src/Debug.cc +++ b/AMDiS/src/Debug.cc @@ -178,6 +178,21 @@ namespace AMDiS { } + Element* getLevel0ParentElement(Mesh *mesh, int elIndex) + { + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); + while (elInfo) { + if (elInfo->getElement()->getIndex() == elIndex) + return elInfo->getMacroElement()->getElement(); + + elInfo = stack.traverseNext(elInfo); + } + + return NULL; + } + + Element* getParentElement(Mesh *mesh, Element *el) { TraverseStack stack; @@ -243,6 +258,22 @@ namespace AMDiS { } + void printElementCoords(const FiniteElemSpace *feSpace, Element *el) + { + FUNCNAME("debug:printElementCoords()"); + + Mesh *mesh = feSpace->getMesh(); + + for (int i = 0; i <= mesh->getDim(); i++) { + DegreeOfFreedom dof = el->getDOF(i, 0); + WorldVector<double> coords; + mesh->getDofIndexCoords(dof, feSpace, coords); + MSG("%d-th DOF of element %d: %f %f %f\n", + i, el->getIndex(), coords[0], coords[1], coords[2]); + } + } + + void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof) { FUNCNAME("debug::printInfoByDof()"); @@ -557,6 +588,48 @@ namespace AMDiS { } + void importDofVectorByCoords(DOFVector<double>* vec, std::string filename) + { + DOFVector<WorldVector<double> > coords(vec->getFeSpace(), "dofCoords"); + vec->getFeSpace()->getMesh()->getDofIndexCoords(vec->getFeSpace(), coords); + int dim = vec->getFeSpace()->getMesh()->getDim(); + + std::ifstream file; + file.open(filename.c_str()); + + int nEntries; + file >> nEntries; + + for (int i = 0; i < nEntries; i++) { + WorldVector<double> dofCoords; + for (int j = 0; j < dim; j++) + file >> dofCoords[j]; + + double value; + file >> value; + + + DOFIterator<WorldVector<double> > it(&coords, USED_DOFS); + for (it.reset(); !it.end(); ++it) { + bool found = true; + for (int j = 0; j < dim; j++) { + if (fabs((*it)[j] - dofCoords[j]) > 1e-8) { + found = false; + break; + } + } + + if (found) { + (*vec)[it.getDOFIndex()] = value; + break; + } + } + } + + file.close(); + } + + void exportDofVectorByCoords(const DOFVector<double>* vec, std::string filename) { diff --git a/AMDiS/src/Debug.h b/AMDiS/src/Debug.h index 3112e75eac549dd8e06d8026f0e61cb092ddaf0e..bdc7cfd4f6bc64bd59d28cea6800b7974f3942e1 100644 --- a/AMDiS/src/Debug.h +++ b/AMDiS/src/Debug.h @@ -86,6 +86,8 @@ namespace AMDiS { Element* getLevel0ParentElement(Mesh *mesh, Element *el); + Element* getLevel0ParentElement(Mesh *mesh, int elIndex); + Element* getParentElement(Mesh *mesh, Element *el); Element* getParentElement(Mesh *mesh, int elIndex); @@ -93,6 +95,8 @@ namespace AMDiS { Element* getElement(Mesh *mesh, int elIndex); void printElementInfo(Element *el); + + void printElementCoords(const FiniteElemSpace *feSpace, Element *el); void printInfoByDof(FiniteElemSpace *feSpace, DegreeOfFreedom dof); @@ -141,6 +145,8 @@ namespace AMDiS { int getLocalNeighbourIndex(Mesh *mesh, int elIndex, int neighIndex); + void importDofVectorByCoords(DOFVector<double>* vec, std::string filename); + void exportDofVectorByCoords(const DOFVector<double>* vec, std::string filename); diff --git a/AMDiS/src/MatrixVector.h b/AMDiS/src/MatrixVector.h index ef70750889b0efb5edff8f1f6b68e77e489de72d..785a54fc2ea0e179f919446340cbd9ba5d952d8c 100644 --- a/AMDiS/src/MatrixVector.h +++ b/AMDiS/src/MatrixVector.h @@ -171,7 +171,7 @@ namespace AMDiS { /// Access to the i-th vector element for const vectors. inline const T& operator[] (int i) const { - TEST_EXIT_DBG(i < size && i >= 0)("invalid index\n"); + TEST_EXIT_DBG(i < size && i >= 0)("Invalid index %d!\n", i); return valArray[i]; } diff --git a/AMDiS/src/Tetrahedron.cc b/AMDiS/src/Tetrahedron.cc index c812172c0cec20f8c4165e55868a89d72484b623..8157c120ee11afc30c1ffb4b4e64ab36771985cb 100644 --- a/AMDiS/src/Tetrahedron.cc +++ b/AMDiS/src/Tetrahedron.cc @@ -214,6 +214,10 @@ namespace AMDiS { BoundaryObject bound, DofContainer& dofs) const { + FUNCNAME("Tetrahedron::getVertexDofsAtFace()"); + + // MSG("go at el %d with face %d\n", this->getIndex(), bound.ithObj); + if (!child[0]) return; @@ -246,7 +250,7 @@ namespace AMDiS { for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); it != bound.excludedSubstructures.end(); ++it) if (it->first == EDGE && it->second == 0) - addDof = false; + addDof = false; if (bound.reverseMode) { child[1]->getVertexDofs(feSpace, nextBound1, dofs); @@ -327,8 +331,8 @@ namespace AMDiS { case FACE: for (ExcludeList::iterator it = bound.excludedSubstructures.begin(); it != bound.excludedSubstructures.end(); ++it) { - if (it->first == EDGE) - it->second = edgeOfChild[bound.elType][ithChild][it->second]; + if (it->first == EDGE && it->second != -1) + it->second = edgeOfChild[bound.elType][ithChild][it->second]; } bound.ithObj = sideOfChild[bound.elType][ithChild][bound.ithObj]; diff --git a/AMDiS/src/ZeroOrderAssembler.cc b/AMDiS/src/ZeroOrderAssembler.cc index 845ca342012bc048321150e93bf79ff1b025176f..3e154f5e5937ac5209f605067897d3030bef39dd 100644 --- a/AMDiS/src/ZeroOrderAssembler.cc +++ b/AMDiS/src/ZeroOrderAssembler.cc @@ -68,9 +68,9 @@ namespace AMDiS { newAssembler = new StandardZOA(op, assembler, quad); } else { if (pwConst) - newAssembler = new PrecalcZOA(op, assembler, quad); + newAssembler = new PrecalcZOA(op, assembler, quad); else - newAssembler = new FastQuadZOA(op, assembler, quad); + newAssembler = new FastQuadZOA(op, assembler, quad); } subAssemblers->push_back(newAssembler); diff --git a/AMDiS/src/ZeroOrderTerm.cc b/AMDiS/src/ZeroOrderTerm.cc index 67e485aa63484afe54b4a613d940d84dc307bccf..94a8b112dc1c74b20981a6b352f5fecee384a178 100644 --- a/AMDiS/src/ZeroOrderTerm.cc +++ b/AMDiS/src/ZeroOrderTerm.cc @@ -36,7 +36,7 @@ namespace AMDiS { { if (f) { for (int iq = 0; iq < nPoints; iq++) - C[iq] += (*f)(vecAtQPs[iq]); + C[iq] += (*f)(vecAtQPs[iq]); } else { for (int iq = 0; iq < nPoints; iq++) C[iq] += vecAtQPs[iq]; diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index 34ceb5fd90add257cc41448e514191bc4912fd1b..50ca44af551810a601d6a69bfcfed7892c0f1fe9 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -147,6 +147,7 @@ namespace AMDiS { debug::testSortedDofs(mesh, elMap); ParallelDebug::testInteriorBoundary(*this); ParallelDebug::testCommonDofs(*this, true); + ParallelDebug::testGlobalIndexByCoords(*this); MSG("Debug mode tests finished!\n"); debug::writeMesh(feSpace, -1, "macro_mesh"); @@ -595,19 +596,6 @@ namespace AMDiS { if (elCode.getCode() != recvCodes[i].getCode()) { TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n"); -#if 0 - MSG("CHECK EL %d %d %d %d WITH El %d %d %d %d in RANK %d\n", - boundIt->rankObj.elIndex, - boundIt->rankObj.subObj, - boundIt->rankObj.ithObj, - boundIt->rankObj.elType, - boundIt->neighObj.elIndex, - boundIt->neighObj.subObj, - boundIt->neighObj.ithObj, - boundIt->neighObj.elType, - it->first); -#endif - bool b = startFitElementToMeshCode(recvCodes[i], boundIt->rankObj.el, boundIt->rankObj.subObj, @@ -699,7 +687,6 @@ namespace AMDiS { refineManager->setMesh(el->getMesh()); refineManager->setStack(&stack); refineManager->refineFunction(elInfo); - MSG("REFINE EL %d %d %d\n", elInfo->getElement()->getIndex(), el->getChild(0)->getIndex(), el->getChild(1)->getIndex()); meshChanged = true; } @@ -778,7 +765,6 @@ namespace AMDiS { refineManager->setMesh(el->getMesh()); refineManager->setStack(&stack); refineManager->refineFunction(elInfo); - MSG("REFINE EL %d %d %d\n", elInfo->getElement()->getIndex(), el->getChild(0)->getIndex(), el->getChild(1)->getIndex()); meshChanged = true; } @@ -1221,9 +1207,14 @@ namespace AMDiS { bound.rankObj.subObj = geoIndex; bound.rankObj.ithObj = rankBoundEl.ithObject; - if (geoIndex == FACE) - for (int edgeNo = 0; edgeNo < 6; edgeNo++) - bound.rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo)); + if (geoIndex == FACE) { + for (int edgeNo = 0; edgeNo < 3; edgeNo++) { + int edgeOfFace = + bound.rankObj.el->getEdgeOfFace(bound.rankObj.ithObj, edgeNo); + + bound.rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeOfFace)); + } + } if (owner == mpiRank) { @@ -1601,7 +1592,7 @@ namespace AMDiS { it->rankObj.el->getNonVertexDofs(feSpace, it->rankObj, dofs); for (unsigned int i = 0; i < dofs.size(); i++) - sendDofs[it.getRank()].push_back(dofs[i]); + sendDofs[it.getRank()].push_back(dofs[i]); } @@ -1638,7 +1629,7 @@ namespace AMDiS { // the dofs in ranks partition that are owned by other rank are set to false. isRankDof.clear(); for (int i = 0; i < nRankAllDofs; i++) - isRankDof[i] = true; + isRankDof[i] = true; // Stores for all rank owned dofs a new global index. DofIndexMap rankDofsNewGlobalIndex; @@ -1697,8 +1688,14 @@ namespace AMDiS { lastMeshChangeIndex = mesh->getChangeIndex(); #if (DEBUG != 0) + std::stringstream oss; + oss << "elementIndex-" << mpiRank << ".vtu"; + debug::writeElementIndexMesh(mesh, oss.str()); + debug::testSortedDofs(mesh, elMap); ParallelDebug::testCommonDofs(*this, true); + ParallelDebug::testGlobalIndexByCoords(*this); + ParallelDebug::writeDebugFile(*this, "mpi-dbg", "dat"); #if 0 MSG("------------- Debug information -------------\n"); @@ -1721,16 +1718,6 @@ namespace AMDiS { MSG("%d recv DOF: %d\n", i, *(it->second[i])); } */ - -#if 0 - std::set<const DegreeOfFreedom*> testDofs; - debug::getAllDofs(feSpace, testDofs); - for (std::set<const DegreeOfFreedom*>::iterator it = testDofs.begin(); - it != testDofs.end(); ++it) { - MSG("DOF %d: mapLocalGlobalDofs = %d isRankDof = %d\n", **it, - mapLocalGlobalDofs[**it], isRankDof[**it]); - } -#endif #endif #endif } diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc index d058d0252989ac77377ec5faf241a4568f86ecf9..15289a95c63327128c673875009857ca927ffaa3 100644 --- a/AMDiS/src/parallel/ParallelDebug.cc +++ b/AMDiS/src/parallel/ParallelDebug.cc @@ -265,6 +265,59 @@ namespace AMDiS { } + void ParallelDebug::testGlobalIndexByCoords(MeshDistributor &pdb) + { + FUNCNAME("ParallelDebug::testGlobalIndexByCoords()"); + + DOFVector<WorldVector<double> > coords(pdb.feSpace, "tmp"); + pdb.mesh->getDofIndexCoords(pdb.feSpace, coords); + + typedef std::map<WorldVector<double>, int> CoordsIndexMap; + CoordsIndexMap coordsToIndex; + + DOFIterator<WorldVector<double> > it(&coords, USED_DOFS); + for (it.reset(); !it.end(); ++it) + coordsToIndex[(*it)] = pdb.mapLocalGlobalDofs[it.getDOFIndex()]; + + StdMpi<CoordsIndexMap> stdMpi(pdb.mpiComm, true); + for (RankToDofContainer::iterator it = pdb.sendDofs.begin(); + it != pdb.sendDofs.end(); ++it) + stdMpi.send(it->first, coordsToIndex); + for (RankToDofContainer::iterator it = pdb.recvDofs.begin(); + it != pdb.recvDofs.end(); ++it) + stdMpi.recv(it->first); + + stdMpi.startCommunication<double>(MPI_DOUBLE); + + int foundError = 0; + for (RankToDofContainer::iterator it = pdb.recvDofs.begin(); + it != pdb.recvDofs.end(); ++it) { + CoordsIndexMap& otherCoords = stdMpi.getRecvData(it->first); + + for (CoordsIndexMap::iterator coordsIt = otherCoords.begin(); + coordsIt != otherCoords.end(); ++coordsIt) { + if (coordsToIndex.count(coordsIt->first) == 1 && + coordsToIndex[coordsIt->first] != coordsIt->second) { + std::stringstream oss; + oss.precision(5); + oss << "DOF at coords "; + for (int i = 0; i < Global::getGeo(WORLD); i++) + oss << coordsIt->first[i] << " "; + oss << " do not fit together on rank " + << pdb.getMpiRank() << " (global index: " + << coordsToIndex[coordsIt->first] << " and on rank " + << it->first << " (global index: " << coordsIt->second << ")"; + + MSG("[DBG] %s\n", oss.str().c_str()); + foundError = 1; + } + } + } + + mpi::globalAdd(foundError); + TEST_EXIT(foundError == 0)("Error found on at least on rank!\n"); + } + void ParallelDebug::testAllElements(MeshDistributor &pdb) { FUNCNAME("ParallelDebug::testAllElements()"); @@ -467,8 +520,8 @@ namespace AMDiS { } - void ParallelDebug::writeCoordsFile(MeshDistributor &pdb, - std::string prefix, std::string postfix) + void ParallelDebug::writeDebugFile(MeshDistributor &pdb, + std::string prefix, std::string postfix) { FUNCNAME("ParallelDebug::writeCoordsFile()"); @@ -478,16 +531,45 @@ namespace AMDiS { DOFVector<WorldVector<double> > coords(pdb.feSpace, "tmp"); pdb.mesh->getDofIndexCoords(pdb.feSpace, coords); + typedef std::map<int, std::vector<DegreeOfFreedom> > ElDofMap; + ElDofMap elDofMap; + TraverseStack stack; + const BasisFunction *basisFcts = pdb.feSpace->getBasisFcts(); + std::vector<DegreeOfFreedom> localIndices(basisFcts->getNumber()); + ElInfo *elInfo = stack.traverseFirst(pdb.mesh, -1, Mesh::CALL_LEAF_EL); + while (elInfo) { + basisFcts->getLocalIndices(elInfo->getElement(), + pdb.feSpace->getAdmin(), localIndices); + elDofMap[elInfo->getElement()->getIndex()] = localIndices; + elInfo = stack.traverseNext(elInfo); + } + + // === Write informations about all DOFs. === + std::ofstream file; file.open(filename.str().c_str()); file << coords.getUsedSize() << "\n"; DOFIterator<WorldVector<double> > it(&coords, USED_DOFS); for (it.reset(); !it.end(); ++it) { - file << pdb.mapLocalGlobalDofs[it.getDOFIndex()]; + file << it.getDOFIndex() << " " + << pdb.mapLocalGlobalDofs[it.getDOFIndex()] << " " + << pdb.getIsRankDof(it.getDOFIndex()); for (int i = 0; i < pdb.mesh->getDim(); i++) file << " " << (*it)[i]; file << "\n"; } + + // === Write to all elements in ranks mesh the included dofs. === + + file << elDofMap.size() << "\n"; + file << basisFcts->getNumber() << "\n"; + for (ElDofMap::iterator it = elDofMap.begin(); it != elDofMap.end(); ++it) { + file << it->first << "\n"; + for (int i = 0; i < basisFcts->getNumber(); i++) + file << it->second[i] << " "; + file << "\n"; + } + file.close(); } diff --git a/AMDiS/src/parallel/ParallelDebug.h b/AMDiS/src/parallel/ParallelDebug.h index 1318f442c5775ffa888e4ea2d118590b77f03c51..86d26467aafa4f7b875306548b7c75f46c66e360 100644 --- a/AMDiS/src/parallel/ParallelDebug.h +++ b/AMDiS/src/parallel/ParallelDebug.h @@ -42,7 +42,7 @@ namespace AMDiS { /** \brief * This function is used for debugging only. It traverses all interior boundaries - * and compares the dof indices on them with the dof indices of the boundarys + * and compares the DOF indices on them with the dof indices of the boundarys * neighbours. The function fails, when dof indices on an interior boundary do * not fit together. * @@ -52,6 +52,16 @@ namespace AMDiS { */ static void testCommonDofs(MeshDistributor &pdb, bool printCoords = false); + /** \brief + * This function is used for debugging only. It checks if on all ranks DOFs with + * the same coordinates have the same global index. For this, the function genartes + * on all ranks a list of all DOFs with their coordinates and global indices and + * sends the list to all neighbouring ranks. + * + * \param[in] pdb Parallel problem definition used for debugging. + */ + static void testGlobalIndexByCoords(MeshDistributor &pdb); + /** \brief * Tests if all elements in the macro mesh are memeber of exactly one rank. * @@ -114,8 +124,8 @@ namespace AMDiS { */ static void printBoundaryInfo(MeshDistributor &pdb); - static void writeCoordsFile(MeshDistributor &pdb, - std::string prefix, std::string postfix); + static void writeDebugFile(MeshDistributor &pdb, + std::string prefix, std::string postfix); }; } // namespace AMDiS diff --git a/AMDiS/src/parallel/PetscSolver.cc b/AMDiS/src/parallel/PetscSolver.cc index 45bed7ed4c2e4a41b56e138b45c922e9bc3e184b..4f2eaaf98786c45d5419bffb74b1055c358a0c04 100644 --- a/AMDiS/src/parallel/PetscSolver.cc +++ b/AMDiS/src/parallel/PetscSolver.cc @@ -30,7 +30,7 @@ namespace AMDiS { clock_t first = clock(); fillPetscMatrix(systemMatrix, rhs); - solvePetscMatrix(*solution, adaptInfo); + solvePetscMatrix(*solution, adaptInfo); #ifdef _OPENMP INFO(info, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n", @@ -94,7 +94,7 @@ namespace AMDiS { // Calculate the exact position of the column index in the petsc matrix. int colIndex = globalColDof * dispMult + dispAddCol; - // Set only non null values. + // Set only non zero values. if (value(*icursor) != 0.0 || rowIndex == colIndex) { // If the current row is not periodic, but the current dof index is periodic, @@ -116,7 +116,7 @@ namespace AMDiS { values.push_back(value(*icursor) * scalFactor); } - } else { + } else { // The col dof index is not periodic, simple add entry. cols.push_back(colIndex); values.push_back(value(*icursor)); @@ -409,7 +409,6 @@ namespace AMDiS { // === Create PETSc matrix with the computed nnz data structure. === MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows, - // 0, PETSC_NULL, 0, PETSC_NULL, &petscMatrix); 0, d_nnz, 0, o_nnz, &petscMatrix); #if (DEBUG != 0) @@ -481,10 +480,6 @@ namespace AMDiS { // Do not delete the solution vector, use it for the initial guess. // KSPSetInitialGuessNonzero(solver, PETSC_TRUE); -#if (DEBUG != 0) - // ParallelDomainDbg::writeCoordsFile(*meshDistributor, "mpi-coords", "dat"); -#endif - // === Run Petsc. === diff --git a/AMDiS/src/parallel/StdMpi.cc b/AMDiS/src/parallel/StdMpi.cc index fd1c953648b514c6f8d714222d22987ff783d3f5..73c93978ddb341b7d49f58d898d8d3a3927eae20 100644 --- a/AMDiS/src/parallel/StdMpi.cc +++ b/AMDiS/src/parallel/StdMpi.cc @@ -12,6 +12,11 @@ namespace AMDiS { return data.size() * 2; } + int intSizeOf(std::map<WorldVector<double>, int> &data) + { + return data.size() * (Global::getGeo(WORLD) + 1); + } + int intSizeOf(std::vector<MeshStructure> &data) { int s = 0; @@ -91,6 +96,38 @@ namespace AMDiS { data.push_back(std::make_pair(buf[i * 2], buf[i * 2 + 1])); } + void makeBuf(std::map<WorldVector<double>, int> &data, double* buf) + { + int i = 0; + for (std::map<WorldVector<double>, int>::iterator it = data.begin(); + it != data.end(); ++it) { + for (int j = 0; j < Global::getGeo(WORLD); j++) + buf[i++] = it->first[j]; + buf[i++] = static_cast<double>(it->second); + } + } + + void makeFromBuf(std::map<WorldVector<double>, int> &data, double* buf, int bufSize) + { + if (bufSize == 0) + return; + + int oneEntrySize = Global::getGeo(WORLD) + 1; + int nEntry = bufSize / oneEntrySize; + + TEST_EXIT(bufSize % oneEntrySize == 0)("This should not happen!\n"); + + data.clear(); + int i = 0; + WorldVector<double> coords; + + for (int j = 0; j < nEntry; j++) { + for (int k = 0; k < Global::getGeo(WORLD); k++) + coords[k] = buf[i++]; + data[coords] = static_cast<int>(buf[i++]); + } + } + void makeBuf(std::vector<MeshStructure> &data, unsigned long int *buf) { int pos = 0; diff --git a/AMDiS/src/parallel/StdMpi.h b/AMDiS/src/parallel/StdMpi.h index a5187cb8a3ab65d681cc96a7365dc7c81dd1c26c..379f2a3add2304eeafac98e28581a9968e6eb7b7 100644 --- a/AMDiS/src/parallel/StdMpi.h +++ b/AMDiS/src/parallel/StdMpi.h @@ -33,6 +33,8 @@ namespace AMDiS { int intSizeOf(std::map<const DegreeOfFreedom*, DegreeOfFreedom> &data); + int intSizeOf(std::map<WorldVector<double>, int> &data); + int intSizeOf(std::vector<MeshStructure> &data); int intSizeOf(std::vector<WorldVector<double> > &data); @@ -58,6 +60,10 @@ namespace AMDiS { void makeFromBuf(std::vector<std::pair<DegreeOfFreedom, DegreeOfFreedom> > &data, int *buf, int bufSize); + void makeBuf(std::map<WorldVector<double>, int> &data, double* buf); + + void makeFromBuf(std::map<WorldVector<double>, int> &data, double* buf, int bufSize); + void makeBuf(std::vector<MeshStructure> &data, unsigned long int *buf); void makeFromBuf(std::vector<MeshStructure> &data,