#include #include "ParallelDomainBase.h" #include "ParMetisPartitioner.h" #include "Mesh.h" #include "Traverse.h" #include "ElInfo.h" #include "Element.h" #include "MacroElement.h" #include "PartitionElementData.h" #include "DOFMatrix.h" #include "DOFVector.h" #include "SystemVector.h" #include "VtkWriter.h" #include "ElementDofIterator.h" #include "ProblemStatBase.h" #include "StandardProblemIteration.h" #include "ElementFileWriter.h" #include "VertexVector.h" #include "StdMpi.h" #include "petscksp.h" namespace AMDiS { PetscErrorCode myKSPMonitor(KSP ksp, PetscInt iter, PetscReal rnorm, void *) { if (iter % 100 == 0 && MPI::COMM_WORLD.Get_rank() == 0) std::cout << " Petsc-Iteration " << iter << ": " << rnorm << std::endl; return 0; } inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2) { return (*dof1 < *dof2); } ParallelDomainBase::ParallelDomainBase(ProblemIterationInterface *iIF, ProblemTimeInterface *tIF, FiniteElemSpace *fe, RefinementManager *refineManager) : iterationIF(iIF), timeIF(tIF), name(iIF->getName()), feSpace(fe), mesh(fe->getMesh()), refinementManager(refineManager), initialPartitionMesh(true), nRankDofs(0), rstart(0), nComponents(1), deserialized(false), lastMeshChangeIndex(0) { FUNCNAME("ParallelDomainBase::ParalleDomainBase()"); TEST_EXIT(mesh->getNumberOfDOFAdmin() == 1) ("Only meshes with one DOFAdmin are supported!\n"); TEST_EXIT(mesh->getDOFAdmin(0).getNumberOfPreDOFs(0) == 0) ("Wrong pre dof number for DOFAdmin!\n"); mpiRank = MPI::COMM_WORLD.Get_rank(); mpiSize = MPI::COMM_WORLD.Get_size(); mpiComm = MPI::COMM_WORLD; partitioner = new ParMetisPartitioner(mesh, &mpiComm); } void ParallelDomainBase::initParallelization(AdaptInfo *adaptInfo) { FUNCNAME("ParallelDomainBase::initParallelization()"); TEST_EXIT(mpiSize > 1) ("Parallelization does not work with only one process!\n"); // If the problem has been already read from a file, we do not need to do anything. if (deserialized) return; // Test, if the mesh is the macro mesh only! Paritioning of the mesh is supported // only for macro meshes, so it will not work yet if the mesh is already refined // in some way. testForMacroMesh(); // create an initial partitioning of the mesh partitioner->createPartitionData(); // set the element weights, which are 1 at the very first begin setElemWeights(adaptInfo); // and now partition the mesh partitionMesh(adaptInfo); #if (DEBUG != 0) ElementIdxToDofs elMap; dbgCreateElementMap(elMap); if (mpiRank == 0) writePartitioningMesh("part.vtu"); #endif // === Create new global and local DOF numbering. === // Set of all DOFs of the rank. std::vector rankDofs; // Number of DOFs in ranks partition that are owned by the rank. nRankDofs = 0; // Number of all DOFs in the macro mesh. int nOverallDOFs = 0; createLocalGlobalNumbering(rankDofs, nRankDofs, nOverallDOFs); // === Create interior boundary information === createInteriorBoundaryInfo(); // === Remove all macro elements that are not part of the rank partition. === removeMacroElements(); #if (DEBUG != 0) dbgTestElementMap(elMap); dbgTestInteriorBoundary(); #endif // === Reset all DOFAdmins of the mesh. === updateDofAdmins(); // === Create periodic dof mapping, if there are periodic boundaries. === createPeriodicMap(); // === Global refinements. === int globalRefinement = 0; GET_PARAMETER(0, mesh->getName() + "->global refinements", "%d", &globalRefinement); if (globalRefinement > 0) { refinementManager->globalRefine(mesh, globalRefinement); #if (DEBUG != 0) elMap.clear(); dbgCreateElementMap(elMap); #endif updateLocalGlobalNumbering(nRankDofs, nOverallDOFs); updateDofAdmins(); #if (DEBUG != 0) dbgTestElementMap(elMap); #endif // === Update periodic mapping, if there are periodic boundaries. === createPeriodicMap(); } lastMeshChangeIndex = mesh->getChangeIndex(); #if (DEBUG != 0) dbgTestCommonDofs(true); #endif nRankRows = nRankDofs * nComponents; nOverallRows = nOverallDOFs * nComponents; } void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo) {} void ParallelDomainBase::updateDofAdmins() { int nAdmins = mesh->getNumberOfDOFAdmin(); for (int i = 0; i < nAdmins; i++) { DOFAdmin& admin = const_cast(mesh->getDOFAdmin(i)); for (int j = 0; j < admin.getSize(); j++) admin.setDOFFree(j, true); for (int j = 0; j < static_cast(mapLocalGlobalDOFs.size()); j++) admin.setDOFFree(j, false); admin.setUsedSize(mapLocalGlobalDOFs.size()); admin.setUsedCount(mapLocalGlobalDOFs.size()); admin.setFirstHole(mapLocalGlobalDOFs.size()); } } void ParallelDomainBase::testForMacroMesh() { FUNCNAME("ParallelDomainBase::testForMacroMesh()"); int nMacroElements = 0; TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { TEST_EXIT(elInfo->getLevel() == 0) ("Mesh is already refined! This does not work with parallelization!\n"); nMacroElements++; elInfo = stack.traverseNext(elInfo); } TEST_EXIT(nMacroElements >= mpiSize) ("The mesh has less macro elements than number of mpi processes!\n"); } void ParallelDomainBase::setDofMatrix(DOFMatrix* mat, int dispMult, int dispAddRow, int dispAddCol) { FUNCNAME("ParallelDomainBase::setDofMatrix()"); TEST_EXIT(mat)("No DOFMatrix!\n"); using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits= mtl::traits; typedef DOFMatrix::base_matrix_type Matrix; 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::vector cols; std::vector values; cols.reserve(300); values.reserve(300); // === Traverse all rows of the dof matrix and insert row wise the values === // === to the petsc matrix. === for (cursor_type cursor = begin(mat->getBaseMatrix()), cend = end(mat->getBaseMatrix()); cursor != cend; ++cursor) { cols.clear(); values.clear(); // Global index of the current row dof. int globalRowDof = mapLocalGlobalDOFs[*cursor]; // Test if the current row dof is a periodic dof. bool periodicRow = (periodicDof.count(globalRowDof) > 0); // === Traverse all non zero entries of the row and produce vector cols === // === with the column indices of all row entries and vector values === // === with the corresponding values. for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { // Set only non null values. if (value(*icursor) != 0.0) { // Global index of the current column index. int globalColDof = mapLocalGlobalDOFs[col(*icursor)]; // Calculate the exact position of the column index in the petsc matrix. int colIndex = globalColDof * dispMult + dispAddCol; // If the current row is not periodic, but the current dof index is periodic, // we have to duplicate the value to the other corresponding periodic columns. if (periodicRow == false && periodicDof.count(globalColDof) > 0) { // The value is assign to n matrix entries, therefore, every entry // has only 1/n value of the original entry. double scalFactor = 1.0 / (periodicDof[globalColDof].size() + 1.0); // Insert original entry. cols.push_back(colIndex); values.push_back(value(*icursor) * scalFactor); // Insert the periodic entries. for (std::set::iterator it = periodicDof[globalColDof].begin(); it != periodicDof[globalColDof].end(); ++it) { cols.push_back(*it * dispMult + dispAddCol); values.push_back(value(*icursor) * scalFactor); } } else { // Neigher row nor column dof index is periodic, simple add entry. cols.push_back(colIndex); values.push_back(value(*icursor)); } } } // === Up to now we have assembled on row. Now, the row must be send to the === // === corresponding rows to the petsc matrix. === // Calculate petsc row index. int rowIndex = globalRowDof * dispMult + dispAddRow; if (periodicRow) { // The row dof is periodic, so send dof to all the corresponding rows. double scalFactor = 1.0 / (periodicDof[globalRowDof].size() + 1.0); int diagIndex = -1; for (int i = 0; i < static_cast(values.size()); i++) { // Change only the non diagonal values in the col. For the diagonal test // we compare the global dof indices of the dof matrix (not of the petsc // matrix!). if ((cols[i] - dispAddCol) / dispMult != globalRowDof) values[i] *= scalFactor; else diagIndex = i; } // Send the main row to the petsc matrix. MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); // Set diagonal element to zero, i.e., the diagonal element of the current // row is not send to the periodic row indices. if (diagIndex != -1) values[diagIndex] = 0.0; // Send the row to all periodic row indices. for (std::set::iterator it = periodicDof[globalRowDof].begin(); it != periodicDof[globalRowDof].end(); ++it) { int perRowIndex = *it * dispMult + dispAddRow; MatSetValues(petscMatrix, 1, &perRowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); } } else { // The row dof is not periodic, simply send the row to the petsc matrix. MatSetValues(petscMatrix, 1, &rowIndex, cols.size(), &(cols[0]), &(values[0]), ADD_VALUES); } } } void ParallelDomainBase::setDofVector(Vec& petscVec, DOFVector* vec, int dispMult, int dispAdd) { // Traverse all used dofs in the dof vector. DOFVector::Iterator dofIt(vec, USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { // Calculate global row index of the dof. int globalRow = mapLocalGlobalDOFs[dofIt.getDOFIndex()]; // Calculate petsc index of the row dof. int index = globalRow * dispMult + dispAdd; if (periodicDof.count(globalRow) > 0) { // The dof index is periodic, so devide the value to all dof entries. double value = *dofIt / (periodicDof[globalRow].size() + 1.0); VecSetValues(petscVec, 1, &index, &value, ADD_VALUES); for (std::set::iterator it = periodicDof[globalRow].begin(); it != periodicDof[globalRow].end(); ++it) { index = *it * dispMult + dispAdd; VecSetValues(petscVec, 1, &index, &value, ADD_VALUES); } } else { // The dof index is not periodic. double value = *dofIt; VecSetValues(petscVec, 1, &index, &value, ADD_VALUES); } } } void ParallelDomainBase::fillPetscMatrix(Matrix *mat, SystemVector *vec) { FUNCNAME("ParallelDomainBase::fillPetscMatrix()"); clock_t first = clock(); VecCreate(PETSC_COMM_WORLD, &petscRhsVec); VecSetSizes(petscRhsVec, nRankRows, nOverallRows); VecSetType(petscRhsVec, VECMPI); VecCreate(PETSC_COMM_WORLD, &petscSolVec); VecSetSizes(petscSolVec, nRankRows, nOverallRows); VecSetType(petscSolVec, VECMPI); VecCreate(PETSC_COMM_WORLD, &petscTmpVec); VecSetSizes(petscTmpVec, nRankRows, nOverallRows); VecSetType(petscTmpVec, VECMPI); using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits= mtl::traits; typedef DOFMatrix::base_matrix_type Matrix; int d_nnz[nRankRows]; int o_nnz[nRankRows]; std::map > > sendMatrixEntry; for (int i = 0; i < nRankRows; i++) { d_nnz[i] = 0; o_nnz[i] = 0; } for (int i = 0; i < nComponents; i++) { for (int j = 0; j < nComponents; j++) { if ((*mat)[i][j]) { Matrix bmat = (*mat)[i][j]->getBaseMatrix(); traits::col::type col(bmat); traits::const_value::type value(bmat); typedef traits::range_generator::type cursor_type; typedef traits::range_generator::type icursor_type; for (cursor_type cursor = begin(bmat), cend = end(bmat); cursor != cend; ++cursor) { int r = mapLocalGlobalDOFs[*cursor] * nComponents + i; if (isRankDof[*cursor]) { r -= rstart * nComponents; #if (DEBUG != 0) if (r < 0 || r >= nRankRows) { std::cout << "ERROR in rank: " << mpiRank << std::endl; std::cout << " Wrong r = " << r << " " << *cursor << " " << mapLocalGlobalDOFs[*cursor] << " " << nRankRows << std::endl; ERROR_EXIT("Should not happen!\n"); } #endif for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { if (value(*icursor) != 0.0) { int c = mapLocalGlobalDOFs[col(*icursor)] * nComponents + j; if (c >= rstart * nComponents && c < rstart * nComponents + nRankRows) d_nnz[r]++; else o_nnz[r]++; } } } else { int sendToRank = -1; for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) { for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) { if (**dofIt == *cursor) { sendToRank = it->first; break; } } if (sendToRank != -1) break; } TEST_EXIT_DBG(sendToRank != -1)("Should not happen!\n"); for (icursor_type icursor = begin(cursor), icend = end(cursor); icursor != icend; ++icursor) { if (value(*icursor) != 0.0) { int c = mapLocalGlobalDOFs[col(*icursor)] * nComponents + j; sendMatrixEntry[sendToRank].push_back(std::make_pair(r, c)); } } } } } } } MPI::Request request[sendDofs.size() + recvDofs.size()]; int requestCounter = 0; std::vector sendBuffers; sendBuffers.reserve(recvDofs.size()); for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) { int nSend = sendMatrixEntry[it->first].size(); request[requestCounter++] = mpiComm.Isend(&nSend, 1, MPI_INT, it->first, 0); if (nSend > 0) { int *sendbuf = new int[nSend * 2]; for (int i = 0; i < nSend; i++) { sendbuf[i * 2] = sendMatrixEntry[it->first][i].first; sendbuf[i * 2 + 1] = sendMatrixEntry[it->first][i].second; } sendBuffers.push_back(sendbuf); } else { sendBuffers.push_back(NULL); } } std::vector recvSize; recvSize.reserve(sendDofs.size()); int i = 0; for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it) request[requestCounter++] = mpiComm.Irecv(&(recvSize[i++]), 1, MPI_INT, it->first, 0); MPI::Request::Waitall(requestCounter, request); requestCounter = 0; i = 0; for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) { int nSend = sendMatrixEntry[it->first].size(); if (nSend > 0) request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSend * 2, MPI_INT, it->first, 0); i++; } std::vector recvBuffers; recvBuffers.reserve(sendDofs.size()); i = 0; for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it) { if (recvSize[i] > 0) { int *recvbuf = new int[recvSize[i] * 2]; recvBuffers.push_back(recvbuf); request[requestCounter++] = mpiComm.Irecv(recvbuf, recvSize[i] * 2, MPI_INT, it->first, 0); } else { recvBuffers.push_back(NULL); } i++; } MPI::Request::Waitall(requestCounter, request); for (int j = 0; j < static_cast(sendBuffers.size()); j++) if (sendBuffers[j]) delete [] sendBuffers[j]; i = 0; for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it) { if (recvSize[i] > 0) { for (int j = 0; j < recvSize[i]; j++) { int r = recvBuffers[i][j * 2]; int c = recvBuffers[i][j * 2 + 1]; r -= rstart * nComponents; TEST_EXIT_DBG(r >= 0 && r < nRankRows)("Should not happen!\n"); if (c < rstart * nComponents || c >= rstart * nComponents + nRankRows) o_nnz[r]++; else d_nnz[r]++; } delete [] recvBuffers[i]; } i++; } MatCreateMPIAIJ(PETSC_COMM_WORLD, nRankRows, nRankRows, nOverallRows, nOverallRows, 0, d_nnz, 0, o_nnz, &petscMatrix); #if (DEBUG != 0) int a, b; MatGetOwnershipRange(petscMatrix, &a, &b); TEST_EXIT(a == rstart * nComponents)("Wrong matrix ownership range!\n"); TEST_EXIT(b == rstart * nComponents + nRankRows)("Wrong matrix ownership range!\n"); #endif using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end; namespace traits= mtl::traits; typedef DOFMatrix::base_matrix_type Matrix; for (int i = 0; i < nComponents; i++) for (int j = 0; j < nComponents; j++) if ((*mat)[i][j]) setDofMatrix((*mat)[i][j], nComponents, i, j); MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY); for (int i = 0; i < nComponents; i++) setDofVector(petscRhsVec, vec->getDOFVector(i), nComponents, i); VecAssemblyBegin(petscRhsVec); VecAssemblyEnd(petscRhsVec); INFO(info, 8)("Fill petsc matrix needed %.5f seconds\n", TIME_USED(first, clock())); } void ParallelDomainBase::solvePetscMatrix(SystemVector &vec) { FUNCNAME("ParallelDomainBase::solvePetscMatrix()"); #if 0 // Set old solution to be initiual guess for petsc solver. for (int i = 0; i < nComponents; i++) setDofVector(petscSolVec, vec->getDOFVector(i), nComponents, i); VecAssemblyBegin(petscSolVec); VecAssemblyEnd(petscSolVec); #endif // === Init Petsc solver. === KSP solver; KSPCreate(PETSC_COMM_WORLD, &solver); KSPSetOperators(solver, petscMatrix, petscMatrix, SAME_NONZERO_PATTERN); KSPSetTolerances(solver, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT); KSPSetType(solver, KSPBCGS); KSPMonitorSet(solver, myKSPMonitor, PETSC_NULL, 0); KSPSetFromOptions(solver); // Do not delete the solution vector, use it for the initial guess. // KSPSetInitialGuessNonzero(solver, PETSC_TRUE); // === Run Petsc. === KSPSolve(solver, petscRhsVec, petscSolVec); // === Transfere values from Petsc's solution vectors to the dof vectors. PetscScalar *vecPointer; VecGetArray(petscSolVec, &vecPointer); for (int i = 0; i < nComponents; i++) { DOFVector *dofvec = vec.getDOFVector(i); for (int j = 0; j < nRankDofs; j++) (*dofvec)[mapLocalToDofIndex[j]] = vecPointer[j * nComponents + i]; } VecRestoreArray(petscSolVec, &vecPointer); // === Synchronize dofs at common dofs, i.e., dofs that correspond to more === // === than one partition. === synchVectors(vec); // === Print information about solution process. === int iterations = 0; KSPGetIterationNumber(solver, &iterations); MSG(" Number of iterations: %d\n", iterations); double norm = 0.0; MatMult(petscMatrix, petscSolVec, petscTmpVec); VecAXPY(petscTmpVec, -1.0, petscRhsVec); VecNorm(petscTmpVec, NORM_2, &norm); MSG(" Residual norm: %e\n", norm); // === Destroy Petsc's variables. === MatDestroy(petscMatrix); VecDestroy(petscRhsVec); VecDestroy(petscSolVec); VecDestroy(petscTmpVec); KSPDestroy(solver); } void ParallelDomainBase::synchVectors(SystemVector &vec) { #if 0 StdMpi, std::vector > stdMpi(mpiComm); stdMpi.prepareCommunication(false); for (RankToDofContainer::iterator sendIt = sendDofs.begin(); sendIt != sendDofs.end(); ++sendIt, i++) { std::vector dofs; int nSendDOFs = sendIt->second.size(); dofs.reserve(nComponents * sendIt->second.size()); for (int j = 0; j < nComponents; j++) { DOFVector *dofvec = vec.getDOFVector(j); for (int k = 0; k < nSendDOFs; k++) dofs.push_back((*dofvec)[*((sendIt->second)[k])]); } stdMpi.send(sendIt->first, dofs); } stdMpi.startCommunication(); #endif #if 1 std::vector sendBuffers(sendDofs.size()); std::vector recvBuffers(recvDofs.size()); MPI::Request request[sendDofs.size() + recvDofs.size()]; int requestCounter = 0; int i = 0; for (RankToDofContainer::iterator sendIt = sendDofs.begin(); sendIt != sendDofs.end(); ++sendIt, i++) { int nSendDOFs = sendIt->second.size(); sendBuffers[i] = new double[nSendDOFs * nComponents]; int counter = 0; for (int j = 0; j < nComponents; j++) { DOFVector *dofvec = vec.getDOFVector(j); for (int k = 0; k < nSendDOFs; k++) sendBuffers[i][counter++] = (*dofvec)[*((sendIt->second)[k])]; } request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDOFs * nComponents, MPI_DOUBLE, sendIt->first, 0); } i = 0; for (RankToDofContainer::iterator recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt, i++) { int nRecvDOFs = recvIt->second.size() * nComponents; recvBuffers[i] = new double[nRecvDOFs]; request[requestCounter++] = mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_DOUBLE, recvIt->first, 0); } MPI::Request::Waitall(requestCounter, request); i = 0; for (RankToDofContainer::iterator recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt, i++) { int nRecvDOFs = recvIt->second.size(); int counter = 0; for (int j = 0; j < nComponents; j++) { DOFVector *dofvec = vec.getDOFVector(j); for (int k = 0; k < nRecvDOFs; k++) (*dofvec)[*(recvIt->second)[k]] = recvBuffers[i][counter++]; } delete [] recvBuffers[i]; } for (int i = 0; i < static_cast(sendBuffers.size()); i++) delete [] sendBuffers[i]; #endif } void ParallelDomainBase::checkMeshChange() { if (mesh->getChangeIndex() == lastMeshChangeIndex) return; std::cout << "MESH HAS BEEN CHANGED!" << std::endl; exit(0); lastMeshChangeIndex = mesh->getChangeIndex(); } void ParallelDomainBase::serialize(std::ostream &out, DofContainer &data) { int vecSize = data.size(); SerUtil::serialize(out, vecSize); for (int i = 0; i < vecSize; i++) { int dofIndex = (*data[i]); SerUtil::serialize(out, dofIndex); } } void ParallelDomainBase::deserialize(std::istream &in, DofContainer &data, std::map &dofMap) { FUNCNAME("ParallelDomainBase::deserialize()"); int vecSize = 0; SerUtil::deserialize(in, vecSize); data.resize(vecSize); for (int i = 0; i < vecSize; i++) { int dofIndex = 0; SerUtil::deserialize(in, dofIndex); TEST_EXIT_DBG(dofMap.count(dofIndex) != 0) ("Dof index could not be deserialized correctly!\n"); data[i] = dofMap[dofIndex]; } } void ParallelDomainBase::serialize(std::ostream &out, RankToDofContainer &data) { int mapSize = data.size(); SerUtil::serialize(out, mapSize); for (RankToDofContainer::iterator it = data.begin(); it != data.end(); ++it) { int rank = it->first; SerUtil::serialize(out, rank); serialize(out, it->second); } } void ParallelDomainBase::deserialize(std::istream &in, RankToDofContainer &data, std::map &dofMap) { int mapSize = 0; SerUtil::deserialize(in, mapSize); for (int i = 0; i < mapSize; i++) { int rank = 0; SerUtil::deserialize(in, rank); deserialize(in, data[rank], dofMap); } } double ParallelDomainBase::setElemWeights(AdaptInfo *adaptInfo) { double localWeightSum = 0.0; int elNum = -1; elemWeights.clear(); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { Element *element = elInfo->getElement(); // get partition data PartitionElementData *partitionData = dynamic_cast (element->getElementData(PARTITION_ED)); if (partitionData && partitionData->getPartitionStatus() == IN) { if (partitionData->getLevel() == 0) elNum = element->getIndex(); TEST_EXIT_DBG(elNum != -1)("invalid element number\n"); if (element->isLeaf()) { elemWeights[elNum] += 1.0; localWeightSum += 1.0; } } elInfo = stack.traverseNext(elInfo); } return localWeightSum; } void ParallelDomainBase::partitionMesh(AdaptInfo *adaptInfo) { if (initialPartitionMesh) { initialPartitionMesh = false; partitioner->fillCoarsePartitionVec(&oldPartitionVec); partitioner->partition(&elemWeights, INITIAL); } else { oldPartitionVec = partitionVec; partitioner->partition(&elemWeights, ADAPTIVE_REPART, 100.0 /*0.000001*/); } partitioner->fillCoarsePartitionVec(&partitionVec); } void ParallelDomainBase::createInteriorBoundaryInfo() { FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()"); int nNeighbours = mesh->getGeo(NEIGH); int dim = mesh->getDim(); GeoIndex subObj = CENTER; switch (dim) { case 2: subObj = EDGE; break; case 3: subObj = FACE; break; default: ERROR_EXIT("What is this?\n"); } // === First, traverse the mesh and search for all elements that have an === // === boundary with an element within another partition. === TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH | Mesh::FILL_BOUND); while (elInfo) { Element *element = elInfo->getElement(); PartitionElementData *partitionData = dynamic_cast(element->getElementData(PARTITION_ED)); // Check, if the element is within rank's partition. if (partitionData->getPartitionStatus() == IN) { for (int i = 0; i < nNeighbours; i++) { if (!elInfo->getNeighbour(i)) continue; PartitionElementData *neighbourPartitionData = dynamic_cast(elInfo->getNeighbour(i)-> getElementData(PARTITION_ED)); if (neighbourPartitionData->getPartitionStatus() == OUT) { // We have found an element that is in rank's partition, but has a // neighbour outside of the rank's partition. // === Create information about the boundary between the two elements. === AtomicBoundary bound; bound.rankObj.el = element; bound.rankObj.elIndex = element->getIndex(); bound.rankObj.elType = elInfo->getType(); bound.rankObj.subObj = subObj; bound.rankObj.ithObj = i; // Do not set a pointer to the element, because if will be deleted from // mesh after partitioning and the pointer would become invalid. bound.neighObj.el = NULL; bound.neighObj.elIndex = elInfo->getNeighbour(i)->getIndex(); bound.neighObj.elType = -1; bound.neighObj.subObj = subObj; bound.neighObj.ithObj = elInfo->getSideOfNeighbour(i); if (dim == 2) { // i == 2 => getSideOfNeighbour(i) == 2 TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2) ("Should not happen!\n"); } // Get rank number of the neighbouring element. int otherElementRank = partitionVec[elInfo->getNeighbour(i)->getIndex()]; // === Add the boundary information object to the corresponding overall === // === boundary. There are three possibilities: if the boundary is a === // === periodic boundary, just add it to \ref periodicBounadry. Here it === // === does not matter which rank is responsible for this boundary. === // === Otherwise, the boundary is added either to \ref myIntBoundary or === // === to \ref otherIntBoundary. It dependes on which rank is respon- === // === sible for this boundary. === if (BoundaryManager::isBoundaryPeriodic(elInfo->getBoundary(subObj, i))) { // We have found an element that is at an interior, periodic boundary. AtomicBoundary& b = periodicBoundary.getNewAtomic(otherElementRank); b = bound; } else { // We have found an element that is at an interior, non-periodic boundary. AtomicBoundary& b = ((mpiRank > otherElementRank) ? myIntBoundary.getNewAtomic(otherElementRank) : otherIntBoundary.getNewAtomic(otherElementRank)); b = bound; } } } } elInfo = stack.traverseNext(elInfo); } // === Once we have this information, we must care about the order of the atomic === // === bounds in the three boundary handling object. Eventually all the bound- === // === aries have to be in the same order on both ranks that share the bounday. === std::vector sendBuffers, recvBuffers; // First we communicate myInt/otherIntBoundary, and than the periodic boundaries. MPI::Request request[max(myIntBoundary.boundary.size() + otherIntBoundary.boundary.size(), periodicBoundary.boundary.size())]; int requestCounter = 0; // === The owner of the boundaries send their atomic boundary order to the === // === neighbouring ranks. === for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin(); rankIt != myIntBoundary.boundary.end(); ++rankIt) { int nSendInt = rankIt->second.size(); int* buffer = new int[nSendInt * 2]; for (int i = 0; i < nSendInt; i++) { buffer[i * 2] = (rankIt->second)[i].rankObj.elIndex; buffer[i * 2 + 1] = (rankIt->second)[i].rankObj.ithObj; } sendBuffers.push_back(buffer); request[requestCounter++] = mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0); } for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); rankIt != otherIntBoundary.boundary.end(); ++rankIt) { int nRecvInt = rankIt->second.size(); int *buffer = new int[nRecvInt * 2]; recvBuffers.push_back(buffer); request[requestCounter++] = mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0); } MPI::Request::Waitall(requestCounter, request); // === The information about all neighbouring boundaries has been received. So === // === the rank tests if its own atomic boundaries are in the same order. If === // === not, the atomic boundaries are swaped to the correct order. === int i = 0; for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); rankIt != otherIntBoundary.boundary.end(); ++rankIt) { // === We have received from rank "rankIt->first" the ordered list of element === // === indices. We now have to sort the corresponding list in this rank to === // === get the same order. === for (int j = 0; j < static_cast(rankIt->second.size()); j++) { // If the expected object is not at place, search for it. if ((rankIt->second)[j].neighObj.elIndex != recvBuffers[i][j * 2] || (rankIt->second)[j].neighObj.ithObj != recvBuffers[i][j * 2 + 1]) { int k = j + 1; for (; k < static_cast(rankIt->second.size()); k++) if ((rankIt->second)[k].neighObj.elIndex == recvBuffers[i][j * 2] && (rankIt->second)[k].neighObj.ithObj == recvBuffers[i][j * 2 + 1]) break; // The element must always be found, because the list is just in another order. TEST_EXIT_DBG(k < static_cast(rankIt->second.size())) ("Should never happen!\n"); // Swap the current with the found element. AtomicBoundary tmpBound = (rankIt->second)[k]; (rankIt->second)[k] = (rankIt->second)[j]; (rankIt->second)[j] = tmpBound; } } delete [] recvBuffers[i]; i++; } for (int i = 0; i < static_cast(sendBuffers.size()); i++) delete [] sendBuffers[i]; recvBuffers.clear(); sendBuffers.clear(); // === Do the same for the periodic boundaries. === requestCounter = 0; for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin(); rankIt != periodicBoundary.boundary.end(); ++rankIt) { TEST_EXIT_DBG(rankIt->first != mpiRank) ("It is no possible to have an interior boundary within a rank partition!\n"); if (rankIt->first < mpiRank) { int nSendInt = rankIt->second.size(); int* buffer = new int[nSendInt * 2]; for (int i = 0; i < nSendInt; i++) { buffer[i * 2] = (rankIt->second)[i].rankObj.elIndex; buffer[i * 2 + 1] = (rankIt->second)[i].rankObj.ithObj; } sendBuffers.push_back(buffer); request[requestCounter++] = mpiComm.Isend(buffer, nSendInt * 2, MPI_INT, rankIt->first, 0); } else { int nRecvInt = rankIt->second.size(); int *buffer = new int[nRecvInt * 2]; recvBuffers.push_back(buffer); request[requestCounter++] = mpiComm.Irecv(buffer, nRecvInt * 2, MPI_INT, rankIt->first, 0); } } MPI::Request::Waitall(requestCounter, request); i = 0; for (RankToBoundMap::iterator rankIt = periodicBoundary.boundary.begin(); rankIt != periodicBoundary.boundary.end(); ++rankIt) { if (rankIt->first > mpiRank) { for (int j = 0; j < static_cast(rankIt->second.size()); j++) if (periodicBoundary.boundary[rankIt->first][j].neighObj.elIndex != recvBuffers[i][j * 2] || periodicBoundary.boundary[rankIt->first][j].neighObj.ithObj != recvBuffers[i][j * 2 + 1]) { int k = j + 1; for (; k < static_cast(rankIt->second.size()); k++) if (periodicBoundary.boundary[rankIt->first][k].neighObj.elIndex == recvBuffers[i][j * 2] && periodicBoundary.boundary[rankIt->first][k].neighObj.ithObj == recvBuffers[i][j * 2 + 1]) break; // The element must always be found, because the list is just in another order. TEST_EXIT_DBG(k < static_cast(rankIt->second.size())) ("Should never happen!\n"); // Swap the current with the found element. AtomicBoundary tmpBound = (rankIt->second)[k]; (rankIt->second)[k] = (rankIt->second)[j]; (rankIt->second)[j] = tmpBound; } delete [] recvBuffers[i]; i++; } } for (int i = 0; i < static_cast(sendBuffers.size()); i++) delete [] sendBuffers[i]; } void ParallelDomainBase::removeMacroElements() { std::vector macrosToRemove; for (std::deque::iterator it = mesh->firstMacroElement(); it != mesh->endOfMacroElements(); ++it) { PartitionElementData *partitionData = dynamic_cast ((*it)->getElement()->getElementData(PARTITION_ED)); if (partitionData->getPartitionStatus() != IN) macrosToRemove.push_back(*it); } mesh->removeMacroElements(macrosToRemove, feSpace); } void ParallelDomainBase::createLocalGlobalNumbering(DofContainer& rankDofs, int& nRankDofs, int& nOverallDOFs) { FUNCNAME("ParallelDomainBase::createLocalGlobalNumbering()"); // === Get rank information about DOFs. === // Stores to each DOF pointer the set of ranks the DOF is part of. std::map > partitionDOFs; DofContainer rankAllDofs; DofToRank boundaryDofs; createDofMemberInfo(partitionDOFs, rankDofs, rankAllDofs, boundaryDofs, vertexDof); nRankDofs = rankDofs.size(); nOverallDOFs = partitionDOFs.size(); // === Get starting position for global rank dof ordering. ==== rstart = 0; mpiComm.Scan(&nRankDofs, &rstart, 1, MPI_INT, MPI_SUM); rstart -= nRankDofs; // === Create for all dofs in rank new indices. The new index must start at === // === index 0, must be continues and have the same order as the indices === // === had before. === // Do not change the indices now, but create a new indexing and store it here. DofIndexMap rankDofsNewLocalIndex; isRankDof.clear(); int i = 0; for (DofContainer::iterator dofIt = rankAllDofs.begin(); dofIt != rankAllDofs.end(); ++dofIt) { rankDofsNewLocalIndex[*dofIt] = i; // First, we set all dofs in ranks partition to be owend by the rank. Later, // the dofs in ranks partition that are owned by other rank are set to false. isRankDof[i] = true; i++; } // === Create for all rank owned dofs a new global indexing. === // Stores for dofs in rank a new global index. DofIndexMap rankDofsNewGlobalIndex; // Stores for all rank owned dofs a continues local index. DofIndexMap rankOwnedDofsNewLocalIndex; i = 0; for (DofContainer::iterator dofIt = rankDofs.begin(); dofIt != rankDofs.end(); ++dofIt) { rankDofsNewGlobalIndex[*dofIt] = i + rstart; rankOwnedDofsNewLocalIndex[*dofIt] = i; i++; } // === Create information which dof indices must be send and which must === // === be received. === // Stores to each rank a map from DOF pointers (which are all owned by the rank // and lie on an interior boundary) to new global DOF indices. std::map sendNewDofs; // Maps to each rank the number of new DOF indices this rank will receive from. // All these DOFs are at an interior boundary on this rank, but are owned by // another rank. std::map recvNewDofs; for (DofToRank::iterator it = boundaryDofs.begin(); it != boundaryDofs.end(); ++it) { if (it->second == mpiRank) { // If the boundary dof is a rank dof, it must be send to other ranks. // Search for all ranks that have this dof too. for (std::set::iterator itRanks = partitionDOFs[it->first].begin(); itRanks != partitionDOFs[it->first].end(); ++itRanks) { if (*itRanks != mpiRank) { TEST_EXIT_DBG(rankDofsNewGlobalIndex.count(it->first) == 1) ("DOF Key not found!\n"); sendNewDofs[*itRanks][it->first] = rankDofsNewGlobalIndex[it->first]; } } } else { // If the boundary dof is not a rank dof, its new dof index (and later // also the dof values) must be received from another rank. if (recvNewDofs.find(it->second) == recvNewDofs.end()) recvNewDofs[it->second] = 1; else recvNewDofs[it->second]++; } } // === Send and receive the dof indices at boundary. === sendDofs.clear(); for (std::map::iterator sendIt = sendNewDofs.begin(); sendIt != sendNewDofs.end(); ++sendIt) for (DofIndexMap::iterator dofIt = sendIt->second.begin(); dofIt != sendIt->second.end(); ++dofIt) sendDofs[sendIt->first].push_back(dofIt->first); StdMpi > > stdMpi(mpiComm); stdMpi.prepareCommunication(false); stdMpi.send(sendNewDofs); for (std::map::iterator recvIt = recvNewDofs.begin(); recvIt != recvNewDofs.end(); ++recvIt, i++) stdMpi.recv(recvIt->first, recvIt->second * 2); stdMpi.startCommunication(); std::map > > &dofMap = stdMpi.getRecvData(); // === Change dof indices at boundary from other ranks. === // Within this small data structure we track which dof index was already changed. // This is used to avoid the following situation: Assume, there are two dof indices // a and b in boundaryDofs. Then we have to change index a to b and b to c. When // the second rule applies, we have to avoid that not the first b, resulted from // changing a to b, is set to c, but the second one. Therefore, after the first // rule was applied, the dof pointer is set to false in this data structure and // is not allowed to be changed anymore. std::map dofChanged; recvDofs.clear(); for (DofToRank::iterator dofIt = boundaryDofs.begin(); dofIt != boundaryDofs.end(); ++dofIt) dofChanged[dofIt->first] = false; for (std::map > >::iterator rankIt = dofMap.begin(); rankIt != dofMap.end(); ++rankIt) { for (std::vector >::iterator recvDofIt = rankIt->second.begin(); recvDofIt != rankIt->second.end(); ++recvDofIt) { DegreeOfFreedom oldDof = recvDofIt->first; DegreeOfFreedom newGlobalDof = recvDofIt->second; bool found = false; // Iterate over all boundary dofs to find the dof which index we have to change. for (DofToRank::iterator dofIt = boundaryDofs.begin(); dofIt != boundaryDofs.end(); ++dofIt) { if (*(dofIt->first) == oldDof && !dofChanged[dofIt->first]) { dofChanged[dofIt->first] = true; recvDofs[rankIt->first].push_back(dofIt->first); rankDofsNewGlobalIndex[dofIt->first] = newGlobalDof; isRankDof[rankDofsNewLocalIndex[dofIt->first]] = false; found = true; break; } } TEST_EXIT_DBG(found)("Should not happen!\n"); } } // === Create now the local to global index and local to dof index mappings. === createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex, rankDofsNewGlobalIndex); } void ParallelDomainBase::updateLocalGlobalNumbering(int& nRankDofs, int& nOverallDOFs) { FUNCNAME("ParallelDomainBase::updateLocalGlobalNumbering()"); typedef std::set DofSet; // === Get all DOFs in ranks partition. === ElementDofIterator elDofIt(feSpace); DofSet rankDofSet; // The vertexDof list must be recreated from the scratch. Otherwise, it is possible // that it maps dofs, that were removed (this is also possible, if the mesh was // refined, e.g., center dofs of an element are not dofs of the children). DofToBool oldVertexDof = vertexDof; vertexDof.clear(); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { Element *element = elInfo->getElement(); elDofIt.reset(element); do { rankDofSet.insert(elDofIt.getDofPtr()); if (elDofIt.getCurrentPos() == 0) vertexDof[elDofIt.getDofPtr()] = true; else vertexDof[elDofIt.getDofPtr()] = false; } while(elDofIt.next()); elInfo = stack.traverseNext(elInfo); } DofContainer rankAllDofs; for (DofSet::iterator dofIt = rankDofSet.begin(); dofIt != rankDofSet.end(); ++dofIt) rankAllDofs.push_back(*dofIt); sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue); DofContainer rankDofs = rankAllDofs; // === Traverse on interior boundaries and move all not ranked owned DOFs from === // === rankDofs to boundaryDOFs. === RankToDofContainer oldSendDofs = sendDofs; RankToDofContainer oldRecvDofs = recvDofs; sendDofs.clear(); recvDofs.clear(); for (RankToDofContainer::iterator it = oldSendDofs.begin(); it != oldSendDofs.end(); ++it) { for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) { if (oldVertexDof[*dofIt]) sendDofs[it->first].push_back(*dofIt); } } for (RankToDofContainer::iterator it = oldRecvDofs.begin(); it != oldRecvDofs.end(); ++it) { for (DofContainer::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) { if (oldVertexDof[*dofIt]) { recvDofs[it->first].push_back(*dofIt); DofContainer::iterator eraseIt = find(rankDofs.begin(), rankDofs.end(), *dofIt); if (eraseIt != rankDofs.end()) rankDofs.erase(eraseIt); } } } for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin(); it != myIntBoundary.boundary.end(); ++it) { DofContainer &dofsToSend = sendDofs[it->first]; for (std::vector::iterator boundIt = it->second.begin(); boundIt != it->second.end(); ++boundIt) { DofContainer dofs; boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, boundIt->rankObj.elType, dofs); boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, boundIt->rankObj.elType, dofs); for (int i = 0; i < static_cast(dofs.size()); i++) dofsToSend.push_back(dofs[i]); } } for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin(); it != otherIntBoundary.boundary.end(); ++it) { DofContainer &dofsToRecv = recvDofs[it->first]; for (std::vector::iterator boundIt = it->second.begin(); boundIt != it->second.end(); ++boundIt) { DofContainer dofs; boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, boundIt->rankObj.elType, dofs); boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, boundIt->rankObj.elType, dofs); if (mesh->getDim() == 2) { for (int i = static_cast(dofs.size()) - 1; i >= 0; i--) { DofContainer::iterator eraseIt = find(rankDofs.begin(), rankDofs.end(), dofs[i]); if (eraseIt != rankDofs.end()) rankDofs.erase(eraseIt); dofsToRecv.push_back(dofs[i]); } } else { for (int i = 0; i < static_cast(dofs.size()); i++) { DofContainer::iterator eraseIt = find(rankDofs.begin(), rankDofs.end(), dofs[i]); if (eraseIt != rankDofs.end()) rankDofs.erase(eraseIt); dofsToRecv.push_back(dofs[i]); } } } } nRankDofs = rankDofs.size(); // === Get starting position for global rank dof ordering. ==== mpiComm.Scan(&nRankDofs, &rstart, 1, MPI_INT, MPI_SUM); rstart -= nRankDofs; // === Calculate number of overall DOFs of all partitions. === mpiComm.Allreduce(&nRankDofs, &nOverallDOFs, 1, MPI_INT, MPI_SUM); // Do not change the indices now, but create a new indexing a store it here. DofIndexMap rankDofsNewLocalIndex; isRankDof.clear(); int i = 0; for (DofContainer::iterator dofIt = rankAllDofs.begin(); dofIt != rankAllDofs.end(); ++dofIt) { rankDofsNewLocalIndex[*dofIt] = i; // First, we set all dofs in ranks partition to be owend by the rank. Later, // the dofs in ranks partition that are owned by other rank are set to false. isRankDof[i] = true; i++; } // Stores for all rank owned dofs a new global index. DofIndexMap rankDofsNewGlobalIndex; // Stores for all rank owned dofs a continues local index. DofIndexMap rankOwnedDofsNewLocalIndex; i = 0; for (DofContainer::iterator dofIt = rankDofs.begin(); dofIt != rankDofs.end(); ++dofIt) { rankDofsNewGlobalIndex[*dofIt] = i + rstart; rankOwnedDofsNewLocalIndex[*dofIt] = i; i++; } // === Send new DOF indices. === std::vector sendBuffers(sendDofs.size()); std::vector recvBuffers(recvDofs.size()); MPI::Request request[sendDofs.size() + recvDofs.size()]; int requestCounter = 0; i = 0; for (RankToDofContainer::iterator sendIt = sendDofs.begin(); sendIt != sendDofs.end(); ++sendIt, i++) { int nSendDofs = sendIt->second.size(); sendBuffers[i] = new int[nSendDofs]; int c = 0; for (DofContainer::iterator dofIt = sendIt->second.begin(); dofIt != sendIt->second.end(); ++dofIt) sendBuffers[i][c++] = rankDofsNewGlobalIndex[*dofIt]; request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0); } i = 0; for (RankToDofContainer::iterator recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt, i++) { int nRecvDofs = recvIt->second.size(); recvBuffers[i] = new int[nRecvDofs]; request[requestCounter++] = mpiComm.Irecv(recvBuffers[i], nRecvDofs, MPI_INT, recvIt->first, 0); } MPI::Request::Waitall(requestCounter, request); for (int j = 0; j < static_cast(sendBuffers.size()); j++) delete [] sendBuffers[j]; i = 0; for (RankToDofContainer::iterator recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt) { int j = 0; for (DofContainer::iterator dofIt = recvIt->second.begin(); dofIt != recvIt->second.end(); ++dofIt) { rankDofsNewGlobalIndex[*dofIt] = recvBuffers[i][j]; isRankDof[rankDofsNewLocalIndex[*dofIt]] = false; j++; } delete [] recvBuffers[i++]; } // === Create now the local to global index and local to dof index mappings. === createLocalMappings(rankDofsNewLocalIndex, rankOwnedDofsNewLocalIndex, rankDofsNewGlobalIndex); } void ParallelDomainBase::createLocalMappings(DofIndexMap &rankDofsNewLocalIndex, DofIndexMap &rankOwnedDofsNewLocalIndex, DofIndexMap &rankDofsNewGlobalIndex) { mapLocalGlobalDOFs.clear(); mapLocalToDofIndex.clear(); for (DofIndexMap::iterator dofIt = rankDofsNewLocalIndex.begin(); dofIt != rankDofsNewLocalIndex.end(); ++dofIt) { DegreeOfFreedom localDof = dofIt->second; DegreeOfFreedom globalDof = rankDofsNewGlobalIndex[dofIt->first]; *const_cast(dofIt->first) = localDof; mapLocalGlobalDOFs[localDof] = globalDof; } for (DofIndexMap::iterator dofIt = rankOwnedDofsNewLocalIndex.begin(); dofIt != rankOwnedDofsNewLocalIndex.end(); ++dofIt) mapLocalToDofIndex[dofIt->second] = *(dofIt->first); } void ParallelDomainBase::createDofMemberInfo(DofToPartitions& partitionDofs, DofContainer& rankOwnedDofs, DofContainer& rankAllDofs, DofToRank& boundaryDofs, DofToBool& vertexDof) { partitionDofs.clear(); rankOwnedDofs.clear(); rankAllDofs.clear(); boundaryDofs.clear(); vertexDof.clear(); // === Determine to each dof the set of partitions the dof belongs to. === ElementDofIterator elDofIt(feSpace); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { Element *element = elInfo->getElement(); elDofIt.reset(element); do { // Determine to each dof the partition(s) it corresponds to. partitionDofs[elDofIt.getDofPtr()].insert(partitionVec[element->getIndex()]); if (elDofIt.getCurrentPos() == 0) vertexDof[elDofIt.getDofPtr()] = true; else vertexDof[elDofIt.getDofPtr()] = false; } while (elDofIt.next()); elInfo = stack.traverseNext(elInfo); } // === Determine the set of ranks dofs and the dofs ownership at the boundary. === // iterate over all DOFs for (DofToPartitions::iterator it = partitionDofs.begin(); it != partitionDofs.end(); ++it) { bool isInRank = false; // iterate over all partition the current DOF is part of. for (std::set::iterator itpart1 = it->second.begin(); itpart1 != it->second.end(); ++itpart1) { if (*itpart1 == mpiRank) { rankAllDofs.push_back(it->first); if (it->second.size() == 1) { rankOwnedDofs.push_back(it->first); } else { // This dof is at the ranks boundary. It is owned by the rank only if // the rank number is the highest of all ranks containing this dof. bool insert = true; int highestRank = mpiRank; for (std::set::iterator itpart2 = it->second.begin(); itpart2 != it->second.end(); ++itpart2) { if (*itpart2 > mpiRank) insert = false; if (*itpart2 > highestRank) highestRank = *itpart2; } if (insert) rankOwnedDofs.push_back(it->first); boundaryDofs[it->first] = highestRank; } isInRank = true; break; } } if (!isInRank) vertexDof.erase(it->first); } sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue); sort(rankOwnedDofs.begin(), rankOwnedDofs.end(), cmpDofsByValue); } void ParallelDomainBase::createPeriodicMap() { FUNCNAME("ParallelDomainBase::createPeriodicMap()"); // Clear all periodic dof mappings calculated before. We do it from scratch. periodicDof.clear(); MPI::Request request[min(static_cast(periodicBoundary.boundary.size() * 2), 4)]; int requestCounter = 0; std::vector sendBuffers, recvBuffers; // === Each rank traverse its periodic boundaries and sends the dof indices === // === to the rank "on the other side" of the periodic boundary. === for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin(); it != periodicBoundary.boundary.end(); ++it) { TEST_EXIT_DBG(it->first != mpiRank) ("Periodic interior boundary within the rank itself is not possible!\n"); DofContainer dofs; // Create dof indices on the boundary. for (std::vector::iterator boundIt = it->second.begin(); boundIt != it->second.end(); ++boundIt) { boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, boundIt->rankObj.elType, dofs, true); boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, boundIt->rankObj.elType, dofs); } // Send the global indices to the rank on the other side. int *sendbuf = new int[dofs.size()]; for (int i = 0; i < static_cast(dofs.size()); i++) sendbuf[i] = mapLocalGlobalDOFs[*(dofs[i])]; request[requestCounter++] = mpiComm.Isend(sendbuf, dofs.size(), MPI_INT, it->first, 0); sendBuffers.push_back(sendbuf); // Receive from this rank the same number of dofs. int *recvbuf = new int[dofs.size()]; request[requestCounter++] = mpiComm.Irecv(recvbuf, dofs.size(), MPI_INT, it->first, 0); recvBuffers.push_back(recvbuf); } MPI::Request::Waitall(requestCounter, request); // === The rank has received the dofs from the rank on the other side of === // === the boundary. Now it can use them to create the mapping between === // === the periodic dofs in this rank and the corresponding periodic === // === dofs from the other ranks. === std::map > dofFromRank; int i = 0; for (RankToBoundMap::iterator it = periodicBoundary.boundary.begin(); it != periodicBoundary.boundary.end(); ++it) { DofContainer dofs; // Create the dofs on the boundary in inverse order. for (std::vector::iterator boundIt = it->second.begin(); boundIt != it->second.end(); ++boundIt) { DofContainer tmpdofs; boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj.ithObj, boundIt->rankObj.elType, tmpdofs); boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj.ithObj, boundIt->rankObj.elType, tmpdofs, true); for (int j = static_cast(tmpdofs.size()) - 1; j >= 0; j--) dofs.push_back(tmpdofs[j]); } // Added the received dofs to the mapping. for (int j = 0; j < static_cast(dofs.size()); j++) { int globalDofIndex = mapLocalGlobalDOFs[*(dofs[j])]; periodicDof[globalDofIndex].insert(recvBuffers[i][j]); dofFromRank[globalDofIndex].insert(it->first); } delete [] sendBuffers[i]; delete [] recvBuffers[i]; i++; } sendBuffers.clear(); recvBuffers.clear(); if (dofFromRank.size() > 0) TEST_EXIT_DBG(mesh->getDim() == 2) ("Periodic boundary corner problem must be generalized to 3d!\n"); requestCounter = 0; for (std::map >::iterator it = dofFromRank.begin(); it != dofFromRank.end(); ++it) { if (it->second.size() == 2) { TEST_EXIT_DBG(periodicDof[it->first].size() == 2)("Missing periodic dof!\n"); int *sendbuf = new int[2]; sendbuf[0] = *(periodicDof[it->first].begin()); sendbuf[1] = *(++(periodicDof[it->first].begin())); request[requestCounter++] = mpiComm.Isend(sendbuf, 2, MPI_INT, *(it->second.begin()), 0); request[requestCounter++] = mpiComm.Isend(sendbuf, 2, MPI_INT, *(++(it->second.begin())), 0); sendBuffers.push_back(sendbuf); int *recvbuf1 = new int[2]; int *recvbuf2 = new int[2]; request[requestCounter++] = mpiComm.Irecv(recvbuf1, 2, MPI_INT, *(it->second.begin()), 0); request[requestCounter++] = mpiComm.Irecv(recvbuf2, 2, MPI_INT, *(++(it->second.begin())), 0); recvBuffers.push_back(recvbuf1); recvBuffers.push_back(recvbuf2); } } MPI::Request::Waitall(requestCounter, request); i = 0; for (std::map >::iterator it = dofFromRank.begin(); it != dofFromRank.end(); ++it) { if (it->second.size() == 2) { for (int k = 0; k < 2; k++) for (int j = 0; j < 2; j++) if (recvBuffers[i + k][j] != it->first) periodicDof[it->first].insert(recvBuffers[i + k][j]); i++; } } } Flag ParallelDomainBase::oneIteration(AdaptInfo *adaptInfo, Flag toDo) { FUNCNAME("ParallelDomainBase::oneIteration()"); Flag flag = buildAndAdapt(adaptInfo, toDo); if (toDo.isSet(SOLVE)) solve(); if (toDo.isSet(SOLVE_RHS)) ERROR_EXIT("Not yet implemented!\n"); if (toDo.isSet(ESTIMATE)) iterationIF->getProblem()->estimate(adaptInfo); return flag; } Flag ParallelDomainBase::buildAndAdapt(AdaptInfo *adaptInfo, Flag toDo) { FUNCNAME("StandardProblemIteration::buildAndAdapt()"); Flag flag = 0, markFlag = 0; ProblemStatBase *problem = iterationIF->getProblem(); if (toDo.isSet(MARK)) markFlag = problem->markElements(adaptInfo); else markFlag = 3; if (toDo.isSet(BUILD)) problem->buildBeforeRefine(adaptInfo, markFlag); // refine if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED)) flag = problem->refineMesh(adaptInfo); if (toDo.isSet(BUILD)) problem->buildBeforeCoarsen(adaptInfo, markFlag); // coarsen if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED)) flag |= problem->coarsenMesh(adaptInfo); checkMeshChange(); if (toDo.isSet(BUILD)) problem->buildAfterCoarsen(adaptInfo, markFlag, true, true); if (toDo.isSet(BUILD_RHS)) problem->buildAfterCoarsen(adaptInfo, markFlag, false, true); return flag; } void ParallelDomainBase::serialize(std::ostream &out) { SerUtil::serialize(out, elemWeights); SerUtil::serialize(out, initialPartitionMesh); SerUtil::serialize(out, partitionVec); SerUtil::serialize(out, oldPartitionVec); SerUtil::serialize(out, nRankDofs); myIntBoundary.serialize(out); otherIntBoundary.serialize(out); serialize(out, sendDofs); serialize(out, recvDofs); SerUtil::serialize(out, mapLocalGlobalDOFs); SerUtil::serialize(out, mapLocalToDofIndex); SerUtil::serialize(out, isRankDof); serialize(out, vertexDof); serialize(out, periodicDof); SerUtil::serialize(out, rstart); SerUtil::serialize(out, nRankRows); SerUtil::serialize(out, nOverallRows); } void ParallelDomainBase::deserialize(std::istream &in) { SerUtil::deserialize(in, elemWeights); SerUtil::deserialize(in, initialPartitionMesh); SerUtil::deserialize(in, partitionVec); SerUtil::deserialize(in, oldPartitionVec); SerUtil::deserialize(in, nRankDofs); // Create two maps: one from from element indices to the corresponding element // pointers, and one map from Dof indices to the corresponding dof pointers. std::map elIndexMap; std::map dofMap; ElementDofIterator elDofIter(feSpace); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { Element *el = elInfo->getElement(); elIndexMap[el->getIndex()] = el; elInfo = stack.traverseNext(elInfo); elDofIter.reset(el); do { dofMap[elDofIter.getDof()] = elDofIter.getDofPtr(); } while(elDofIter.next()); } myIntBoundary.deserialize(in, elIndexMap); otherIntBoundary.deserialize(in, elIndexMap); deserialize(in, sendDofs, dofMap); deserialize(in, recvDofs, dofMap); SerUtil::deserialize(in, mapLocalGlobalDOFs); SerUtil::deserialize(in, mapLocalToDofIndex); SerUtil::deserialize(in, isRankDof); deserialize(in, vertexDof, dofMap); SerUtil::deserialize(in, rstart); SerUtil::deserialize(in, nRankRows); SerUtil::deserialize(in, nOverallRows); deserialized = true; } void ParallelDomainBase::serialize(std::ostream &out, PeriodicDofMap &data) { int mapSize = data.size(); SerUtil::serialize(out, mapSize); for (PeriodicDofMap::iterator it = data.begin(); it != data.end(); ++it) { DegreeOfFreedom dof = it->first; std::set dofSet = it->second; SerUtil::serialize(out, dof); SerUtil::serialize(out, dofSet); } } void ParallelDomainBase::deserialize(std::istream &in, PeriodicDofMap &data) { data.clear(); int mapSize = 0; SerUtil::deserialize(in, mapSize); for (int i = 0; i < mapSize; i++) { DegreeOfFreedom dof = 0; std::set dofSet; SerUtil::deserialize(in, dof); SerUtil::deserialize(in, dofSet); data[dof] = dofSet; } } void ParallelDomainBase::dbgCreateElementMap(ElementIdxToDofs &elMap) { FUNCNAME("ParallelDomainBase::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: orderDofs(el->getDOF(0), el->getDOF(1), el->getDOF(2), elMap[el->getIndex()]); break; case 3: orderDofs(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 ParallelDomainBase::dbgTestElementMap(ElementIdxToDofs &elMap) { FUNCNAME("ParallelDomainbase::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: orderDofs(el->getDOF(0), el->getDOF(1), el->getDOF(2), vec); break; case 3: orderDofs(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]) { std::cout << "[DBG " << mpiRank << "]: Wrong new dof numeration in element = " << el->getIndex() << std::endl; std::cout << "[DBG " << mpiRank << "]: 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 " << mpiRank << "]: 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 ParallelDomainBase::dbgTestInteriorBoundary() { FUNCNAME("ParallelDomainBase::dbgTestInteriorBoundary()"); std::vector sendBuffers, recvBuffers; MPI::Request request[myIntBoundary.boundary.size() + otherIntBoundary.boundary.size()]; int requestCounter = 0; for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin(); rankIt != myIntBoundary.boundary.end(); ++rankIt) { int nSendInt = rankIt->second.size(); int* buffer = new int[nSendInt]; for (int i = 0; i < nSendInt; i++) buffer[i] = (rankIt->second)[i].rankObj.elIndex; sendBuffers.push_back(buffer); request[requestCounter++] = mpiComm.Isend(buffer, nSendInt, MPI_INT, rankIt->first, 0); } for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); rankIt != otherIntBoundary.boundary.end(); ++rankIt) { int nRecvInt = rankIt->second.size(); int *buffer = new int[nRecvInt]; recvBuffers.push_back(buffer); request[requestCounter++] = mpiComm.Irecv(buffer, nRecvInt, MPI_INT, rankIt->first, 0); } MPI::Request::Waitall(requestCounter, request); for (int i = 0; i < static_cast(sendBuffers.size()); i++) delete [] sendBuffers[i]; int bufCounter = 0; for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); rankIt != otherIntBoundary.boundary.end(); ++rankIt) { TEST_EXIT(rankIt->second.size() == otherIntBoundary.boundary[rankIt->first].size()) ("Boundaries does not fit together!\n"); for (int i = 0; i < static_cast(rankIt->second.size()); i++) { int elIndex1 = recvBuffers[bufCounter][i]; int elIndex2 = otherIntBoundary.boundary[rankIt->first][i].neighObj.elIndex; TEST_EXIT(elIndex1 == elIndex2)("Wrong element index at interior boundary!\n"); } delete [] recvBuffers[bufCounter++]; } } void ParallelDomainBase::dbgTestCommonDofs(bool printCoords) { FUNCNAME("ParallelDomainBase::dbgTestCommonDofs()"); // Maps to each neighbour rank an array of WorldVectors. This array contains the // coordinates of all dofs this rank shares on the interior boundary with the // neighbour rank. A rank sends the coordinates to another rank, if it owns the // boundarys dof. RankToCoords sendCoords; // A rank receives all boundary dofs that are at its interior boundaries but are // not owned by the rank. This map stores for each rank the the coordinates of dofs // this rank expectes to receive from. RankToCoords recvCoords; WorldVector coords; for (RankToDofContainer::iterator it = sendDofs.begin(); it != sendDofs.end(); ++it) { for (std::vector::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) { bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords); if (!b) { std::cout << "ERROR (send): Cannot find DOF in mesh!" << std::endl; std::cout << " mpi rank = " << mpiRank << " send to rank = " << it->first << " dof = " << *dofIt << " = " << **dofIt << std::endl; ERROR_EXIT("Should not happen!\n"); } sendCoords[it->first].push_back(coords); } } for (RankToDofContainer::iterator it = recvDofs.begin(); it != recvDofs.end(); ++it) { for (std::vector::iterator dofIt = it->second.begin(); dofIt != it->second.end(); ++dofIt) { bool b = mesh->getDofIndexCoords(*dofIt, feSpace, coords); if (!b) { std::cout << "ERROR (recv): Cannot find DOF in mesh!" << std::endl; std::cout << " mpi rank = " << mpiRank << " send to rank = " << it->first << " dof = " << *dofIt << " = " << **dofIt << std::endl; ERROR_EXIT("Should not happen!\n"); } recvCoords[it->first].push_back(coords); } } std::vector sendSize(mpiSize, 0); std::vector recvSize(mpiSize, 0); std::vector recvSizeBuffer(mpiSize, 0); MPI::Request request[(mpiSize - 1) * 2]; int requestCounter = 0; for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it) sendSize[it->first] = it->second.size(); for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) recvSize[it->first] = it->second.size(); for (int i = 0; i < mpiSize; i++) { if (i == mpiRank) continue; request[requestCounter++] = mpiComm.Isend(&(sendSize[i]), 1, MPI_INT, i, 0); } for (int i = 0; i < mpiSize; i++) { if (i == mpiRank) continue; request[requestCounter++] = mpiComm.Irecv(&(recvSizeBuffer[i]), 1, MPI_INT, i, 0); } MPI::Request::Waitall(requestCounter, request); for (int i = 0; i < mpiSize; i++) { if (i == mpiRank) continue; if (recvSize[i] != recvSizeBuffer[i]) { std::cout << "Error: MPI rank " << mpiRank << " expectes to receive " << recvSize[i] << " DOFs from rank " << i << ". But this rank sends " << recvSizeBuffer[i] << " DOFs!" << std::endl; ERROR_EXIT("Should not happen!\n"); } } std::map sendCoordsBuffer, recvCoordsBuffer; requestCounter = 0; int dimOfWorld = Global::getGeo(WORLD); for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it) { sendCoordsBuffer[it->first] = new double[it->second.size() * dimOfWorld]; for (int i = 0; i < static_cast(it->second.size()); i++) for (int j = 0; j < dimOfWorld; j++) sendCoordsBuffer[it->first][i * dimOfWorld + j] = (it->second)[i][j]; request[requestCounter++] = mpiComm.Isend(sendCoordsBuffer[it->first], it->second.size() * dimOfWorld, MPI_DOUBLE, it->first, 0); } for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) { recvCoordsBuffer[it->first] = new double[it->second.size() * dimOfWorld]; request[requestCounter++] = mpiComm.Irecv(recvCoordsBuffer[it->first], it->second.size() * dimOfWorld, MPI_DOUBLE, it->first, 0); } MPI::Request::Waitall(requestCounter, request); for (RankToCoords::iterator it = sendCoords.begin(); it != sendCoords.end(); ++it) delete [] sendCoordsBuffer[it->first]; double eps = 1e-13; for (RankToCoords::iterator it = recvCoords.begin(); it != recvCoords.end(); ++it) { for (int i = 0; i < static_cast(it->second.size()); i++) { for (int j = 0; j < dimOfWorld; j++) { bool c = fabs((it->second)[i][j] - recvCoordsBuffer[it->first][i * dimOfWorld + j]) < eps; if (printCoords && !c) { std::cout << "[DBG] i = " << i << std::endl; std::cout.precision(5); std::cout << "[DBG] " << "Rank " << mpiRank << " from rank " << it->first << " expect coords ("; for (int k = 0; k < dimOfWorld; k++) { std::cout << (it->second)[i][k]; if (k + 1 < dimOfWorld) std::cout << " / "; } std::cout << ") received coords ("; for (int k = 0; k < dimOfWorld; k++) { std::cout << recvCoordsBuffer[it->first][i * dimOfWorld + k]; if (k + 1 < dimOfWorld) std::cout << " / "; } std::cout << ")" << std::endl; } TEST_EXIT(c)("Wrong DOFs in rank %d!\n", mpiRank); } } delete [] recvCoordsBuffer[it->first]; } } void ParallelDomainBase::printMapLocalGlobal(int rank) { if (rank == -1 || mpiRank == rank) { std::cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << std::endl; for (DofMapping::iterator it = mapLocalGlobalDOFs.begin(); it != mapLocalGlobalDOFs.end(); it++) { DegreeOfFreedom localdof = -1; if (mapLocalToDofIndex.count(it->first) > 0) localdof = mapLocalToDofIndex[it->first]; std::cout << "DOF " << it->first << " " << it->second << " " << localdof << std::endl; WorldVector coords; mesh->getDofIndexCoords(it->first, feSpace, coords); coords.print(); for (RankToDofContainer::iterator rankit = sendDofs.begin(); rankit != sendDofs.end(); ++rankit) { for (DofContainer::iterator dofit = rankit->second.begin(); dofit != rankit->second.end(); ++dofit) if (**dofit == it->first) std::cout << "SEND DOF TO " << rankit->first << std::endl; } for (RankToDofContainer::iterator rankit = recvDofs.begin(); rankit != recvDofs.end(); ++rankit) { for (DofContainer::iterator dofit = rankit->second.begin(); dofit != rankit->second.end(); ++dofit) if (**dofit == it->first) std::cout << "RECV DOF FROM " << rankit->first << std::endl; } std::cout << "------" << std::endl; } } } void ParallelDomainBase::printMapPeriodic(int rank) { FUNCNAME("ParallelDomainBase::printMapPeriodic()"); if (rank == -1 || mpiRank == rank) { std::cout << "====== DOF MAP PERIODIC ====== " << std::endl; for (PeriodicDofMap::iterator it = periodicDof.begin(); it != periodicDof.end(); ++it) { std::cout << "DOF MAP " << it->first << ": "; for (std::set::iterator dofit = it->second.begin(); dofit != it->second.end(); ++dofit) std::cout << *dofit << " "; std::cout << std::endl; DegreeOfFreedom localdof = -1; for (DofMapping::iterator dofIt = mapLocalGlobalDOFs.begin(); dofIt != mapLocalGlobalDOFs.end(); ++dofIt) if (dofIt->second == it->first) localdof = dofIt->first; TEST_EXIT(localdof != -1)("There is something wrong!\n"); WorldVector coords; mesh->getDofIndexCoords(localdof, feSpace, coords); coords.print(); } } } void ParallelDomainBase::printRankDofs(int rank, DofContainer& rankDofs, DofContainer& rankAllDofs) { if (rank == -1 || mpiRank == rank) { std::cout << "====== RANK DOF INFORMATION ====== " << std::endl; std::cout << " RANK OWNED DOFS: " << std::endl; for (DofContainer::iterator dofit = rankDofs.begin(); dofit != rankDofs.end(); ++dofit) { std::cout << " " << **dofit << std::endl; WorldVector coords; mesh->getDofIndexCoords(*dofit, feSpace, coords); coords.print(); } std::cout << " RANK ALL DOFS: " << std::endl; for (DofContainer::iterator dofit = rankAllDofs.begin(); dofit != rankAllDofs.end(); ++dofit) { std::cout << " " << **dofit << std::endl; WorldVector coords; mesh->getDofIndexCoords(*dofit, feSpace, coords); coords.print(); } } } void ParallelDomainBase::writePartitioningMesh(std::string filename) { FUNCNAME("ParallelDomainBase::writePartitioningMesh()"); std::map vec; TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); while (elInfo) { int index = elInfo->getElement()->getIndex(); vec[index] = partitionVec[index]; elInfo = stack.traverseNext(elInfo); } ElementFileWriter::writeFile(vec, feSpace, filename); } }