#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 "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 << " Iteration " << iter << ": " << rnorm << std::endl; return 0; } inline bool cmpDofsByValue(const DegreeOfFreedom* dof1, const DegreeOfFreedom* dof2) { return (*dof1 < *dof2); } ParallelDomainBase::ParallelDomainBase(const std::string& name, ProblemIterationInterface *iIF, ProblemTimeInterface *tIF, FiniteElemSpace *fe, RefinementManager *refineManager) : iterationIF(iIF), timeIF(tIF), feSpace(fe), mesh(fe->getMesh()), refinementManager(refineManager), initialPartitionMesh(true), nRankDOFs(0), rstart(0), nComponents(1) { 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) { if (mpiSize <= 1) return; #if 0 if (mpiRank == 0) { std::map vec; TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH); while (elInfo) { vec[elInfo->getElement()->getIndex()] = static_cast(elInfo->getElement()->getIndex()); elInfo = stack.traverseNext(elInfo); } ElementFileWriter::writeFile(vec, feSpace, "test.vtu"); } #endif // 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); #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(rankDOFs); // === 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(); #if 0 if (mpiRank == 0) { TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { if (elInfo->getElement()->getIndex() == 4) { WorldVector x; mesh->getDofIndexCoords(elInfo->getElement()->getDOF(0), feSpace, x); std::cout << "FOUND!" << std::endl; x.print(); } elInfo = stack.traverseNext(elInfo); } } if (mpiRank == 1) { TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { if (elInfo->getElement()->getIndex() == 3) { WorldVector x; mesh->getDofIndexCoords(elInfo->getElement()->getDOF(2), feSpace, x); std::cout << "FOUND!" << std::endl; x.print(); } elInfo = stack.traverseNext(elInfo); } } exit(0); #endif // === Global refinements. === int globalRefinement = 0; GET_PARAMETER(0, "testgr", "%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 } #if (DEBUG != 0) DbgTestCommonDofs(true); #endif // === Create petsc matrix. === int nRankRows = nRankDOFs * nComponents; int nOverallRows = nOverallDOFs * nComponents; MatCreate(PETSC_COMM_WORLD, &petscMatrix); MatSetSizes(petscMatrix, nRankRows, nRankRows, nOverallRows, nOverallRows); MatSetType(petscMatrix, MATAIJ); VecCreate(PETSC_COMM_WORLD, &petscRhsVec); VecSetSizes(petscRhsVec, nRankRows, nOverallRows); VecSetType(petscRhsVec, VECMPI); VecCreate(PETSC_COMM_WORLD, &petscSolVec); VecSetSizes(petscSolVec, nRankRows, nOverallRows); VecSetType(petscSolVec, VECMPI); } void ParallelDomainBase::exitParallelization(AdaptInfo *adaptInfo) { MatDestroy(petscMatrix); VecDestroy(petscRhsVec); VecDestroy(petscSolVec); } 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 nor 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::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; 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) if (value(*icursor) != 0.0) { int r = mapLocalGlobalDOFs[row(*icursor)] * dispMult + dispAddRow; int c = mapLocalGlobalDOFs[col(*icursor)] * dispMult + dispAddCol; double v = value(*icursor); MatSetValues(petscMatrix, 1, &r, 1, &c, &v, ADD_VALUES); } } void ParallelDomainBase::setDofVector(DOFVector* vec, int dispMult, int dispAdd) { DOFVector::Iterator dofIt(vec, USED_DOFS); for (dofIt.reset(); !dofIt.end(); ++dofIt) { int index = mapLocalGlobalDOFs[dofIt.getDOFIndex()] * dispMult + dispAdd; double value = *dofIt; VecSetValues(petscRhsVec, 1, &index, &value, ADD_VALUES); } } void ParallelDomainBase::fillPetscMatrix(DOFMatrix *mat, DOFVector *vec) { setDofMatrix(mat); MatAssemblyBegin(petscMatrix, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(petscMatrix, MAT_FINAL_ASSEMBLY); setDofVector(vec); } void ParallelDomainBase::fillPetscMatrix(Matrix *mat, SystemVector *vec) { 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(vec->getDOFVector(i), nComponents, i); } void ParallelDomainBase::solvePetscMatrix(DOFVector *vec) { FUNCNAME("ParallelDomainBase::solvePetscMatrix()"); KSP ksp; PC pc; KSPCreate(PETSC_COMM_WORLD, &ksp); KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN); KSPGetPC(ksp, &pc); // PCSetType(pc, PCNONE); PCSetType(pc, PCILU); KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT); KSPSetType(ksp, KSPBCGS); //KSPSetType(ksp, KSPCG); KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0); KSPSolve(ksp, petscRhsVec, petscSolVec); #if (DEBUG != 0) int size = 0; VecGetLocalSize(petscSolVec, &size); TEST_EXIT(size == nRankDOFs)("Vector and rank DOFs does not fit together!\n"); #endif PetscScalar *vecPointer; VecGetArray(petscSolVec, &vecPointer); for (int i = 0; i < nRankDOFs; i++) (*vec)[mapLocalToDofIndex[i]] = vecPointer[i]; VecRestoreArray(petscSolVec, &vecPointer); 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]; for (int j = 0; j < nSendDOFs; j++) sendBuffers[i][j] = (*vec)[*((sendIt->second)[j])]; request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDOFs, MPI_DOUBLE, sendIt->first, 0); } i = 0; for (RankToDofContainer::iterator recvIt = recvDofs.begin(); recvIt != recvDofs.end(); ++recvIt, i++) { int nRecvDOFs = recvIt->second.size(); 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++) { for (int j = 0; j < static_cast(recvIt->second.size()); j++) (*vec)[*(recvIt->second)[j]] = recvBuffers[i][j]; delete [] recvBuffers[i]; } for (int i = 0; i < static_cast(sendBuffers.size()); i++) delete [] sendBuffers[i]; } void ParallelDomainBase::solvePetscMatrix(SystemVector *vec) { FUNCNAME("ParallelDomainBase::solvePetscMatrix()"); KSP ksp; PC pc; KSPCreate(PETSC_COMM_WORLD, &ksp); KSPSetOperators(ksp, petscMatrix, petscMatrix, DIFFERENT_NONZERO_PATTERN); KSPGetPC(ksp, &pc); // PCSetType(pc, PCNONE); PCSetType(pc, PCJACOBI); KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT); KSPSetType(ksp, KSPBCGS); //KSPSetType(ksp, KSPCG); KSPMonitorSet(ksp, myKSPMonitor, PETSC_NULL, 0); KSPSolve(ksp, petscRhsVec, petscSolVec); #if (DEBUG != 0) int size = 0; VecGetLocalSize(petscSolVec, &size); TEST_EXIT(size == nRankDOFs * nComponents) ("Vector and rank DOFs does not fit together!\n"); #endif 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); 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]; } 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(DofContainer& rankDOFs) { FUNCNAME("ParallelDomainBase::createInteriorBoundaryInfo()"); // === First, create all the information about the interior boundaries. === TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH); while (elInfo) { Element *element = elInfo->getElement(); PartitionElementData *partitionData = dynamic_cast(element->getElementData(PARTITION_ED)); if (partitionData->getPartitionStatus() == IN) { for (int i = 0; i < 3; 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 at an interior boundary. // === Find out, if the boundary part of the element corresponds to the === // === rank or to the rank "on the other side" of the interoir boundary. === const DegreeOfFreedom* boundDOF1 = NULL; const DegreeOfFreedom* boundDOF2 = NULL; switch (i) { case 0: boundDOF1 = element->getDOF(1); boundDOF2 = element->getDOF(2); break; case 1: boundDOF1 = element->getDOF(0); boundDOF2 = element->getDOF(2); break; case 2: boundDOF1 = element->getDOF(0); boundDOF2 = element->getDOF(1); break; default: ERROR_EXIT("Should never happen!\n"); } bool isRankDOF1 = (find(rankDOFs.begin(), rankDOFs.end(), boundDOF1) != rankDOFs.end()); bool isRankDOF2 = (find(rankDOFs.begin(), rankDOFs.end(), boundDOF2) != rankDOFs.end()); bool ranksBoundary = isRankDOF1 || isRankDOF2; #if 0 if (mpiRank == 3 && ranksBoundary && partitionVec[elInfo->getNeighbour(i)->getIndex()] == 2) { std::cout << "ADD MY BOUND " << element->getIndex() << "/" << i << " with " << elInfo->getNeighbour(i)->getIndex() << "/" << elInfo->getSideOfNeighbour(i) << std::endl; } if (mpiRank == 2 && !ranksBoundary && partitionVec[elInfo->getNeighbour(i)->getIndex()] == 3) { std::cout << "ADD OT BOUND " << element->getIndex() << "/" << i << " with " << elInfo->getNeighbour(i)->getIndex() << "/" << elInfo->getSideOfNeighbour(i) << std::endl; } #endif // === And add the part of the interior boundary. === AtomicBoundary& bound = (ranksBoundary ? myIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()]) : otherIntBoundary.getNewAtomicBoundary(partitionVec[elInfo->getNeighbour(i)->getIndex()])); bound.rankObject.el = element; bound.rankObject.elIndex = element->getIndex(); bound.rankObject.subObjAtBoundary = EDGE; bound.rankObject.ithObjAtBoundary = 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.neighbourObject.el = NULL; bound.neighbourObject.elIndex = elInfo->getNeighbour(i)->getIndex(); bound.neighbourObject.subObjAtBoundary = EDGE; bound.neighbourObject.ithObjAtBoundary = elInfo->getSideOfNeighbour(i); // i == 2 => getSideOfNeighbour(i) == 2 TEST_EXIT_DBG(i != 2 || elInfo->getSideOfNeighbour(i) == 2) ("Should not happen!\n"); } } } elInfo = stack.traverseNext(elInfo); } // === Once we have this information, we must care about their order. Eventually === // === all the boundaries have to be in the same order on both ranks that share === // === the bounday. === 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 * 2]; for (int i = 0; i < nSendInt; i++) { buffer[i * 2] = (rankIt->second)[i].rankObject.elIndex; buffer[i * 2 + 1] = (rankIt->second)[i].rankObject.ithObjAtBoundary; } 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); 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].neighbourObject.elIndex != recvBuffers[i][j * 2] || (rankIt->second)[j].neighbourObject.ithObjAtBoundary != recvBuffers[i][j * 2 + 1]) { int k = j + 1; for (; k < static_cast(rankIt->second.size()); k++) if ((rankIt->second)[k].neighbourObject.elIndex == recvBuffers[i][j * 2] && (rankIt->second)[k].neighbourObject.ithObjAtBoundary == 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++]; } 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); 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 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++; } // === 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. === std::vector sendBuffers(sendNewDofs.size()), recvBuffers(recvNewDofs.size()); sendDofs.clear(); recvDofs.clear(); MPI::Request request[sendNewDofs.size() + recvNewDofs.size()]; int requestCounter = 0; i = 0; for (std::map >::iterator sendIt = sendNewDofs.begin(); sendIt != sendNewDofs.end(); ++sendIt, i++) { int nSendDofs = sendIt->second.size() * 2; sendBuffers[i] = new int[nSendDofs]; int c = 0; for (std::map::iterator dofIt = sendIt->second.begin(); dofIt != sendIt->second.end(); ++dofIt) { sendBuffers[i][c++] = *(dofIt->first); sendBuffers[i][c++] = dofIt->second; #if 0 if (mpiRank == 3 && sendIt->first == 2) std::cout << "SEND DOF: " << dofIt->first << std::endl; #endif sendDofs[sendIt->first].push_back(dofIt->first); } request[requestCounter++] = mpiComm.Isend(sendBuffers[i], nSendDofs, MPI_INT, sendIt->first, 0); } i = 0; for (std::map::iterator recvIt = recvNewDofs.begin(); recvIt != recvNewDofs.end(); ++recvIt, i++) { int nRecvDOFs = recvIt->second * 2; recvBuffers[i] = new int[nRecvDOFs]; request[requestCounter++] = mpiComm.Irecv(recvBuffers[i], nRecvDOFs, MPI_INT, recvIt->first, 0); } MPI::Request::Waitall(requestCounter, request); // === Delete send buffers. === for (int j = 0; j < static_cast(sendBuffers.size()); j++) delete [] sendBuffers[j]; // === 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; for (DofToRank::iterator dofIt = boundaryDofs.begin(); dofIt != boundaryDofs.end(); ++dofIt) dofChanged[dofIt->first] = false; i = 0; for (std::map::iterator recvIt = recvNewDofs.begin(); recvIt != recvNewDofs.end(); ++recvIt, i++) { for (int j = 0; j < recvIt->second; j++) { DegreeOfFreedom oldDof = recvBuffers[i][j * 2]; DegreeOfFreedom newGlobalDof = recvBuffers[i][j * 2 + 1]; 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; #if 0 if (mpiRank == 2 && recvIt->first == 3) std::cout << "RECV DOF: " << dofIt->first << std::endl; #endif recvDofs[recvIt->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"); } delete [] recvBuffers[i]; } // === 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; 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()); } while(elDofIt.next()); elInfo = stack.traverseNext(elInfo); } DofContainer rankAllDofs; for (DofSet::iterator dofIt = rankDOFSet.begin(); dofIt != rankDOFSet.end(); ++dofIt) { /* if (mpiRank == 1) { std::cout << "COORDs of dof = " << **dofIt << std::endl; WorldVector x; mesh->getDofIndexCoords(*dofIt, feSpace, x); x.print(); }*/ 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. === sendDofs.clear(); recvDofs.clear(); for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin(); it != myIntBoundary.boundary.end(); ++it) { for (std::vector::iterator boundIt = it->second.begin(); boundIt != it->second.end(); ++boundIt) { #if 0 if (mpiRank == 3 && it->first == 2) std::cout << "GO ON MY BOUND: " << boundIt->rankObject.elIndex << "/" << boundIt->rankObject.ithObjAtBoundary << " with " << boundIt->neighbourObject.elIndex << "/" << boundIt->neighbourObject.ithObjAtBoundary << std::endl; #endif DofContainer dofs; DofContainer &dofsToSend = sendDofs[it->first]; switch (boundIt->rankObject.ithObjAtBoundary) { case 0: dofs.push_back(boundIt->rankObject.el->getDOF(1)); dofs.push_back(boundIt->rankObject.el->getDOF(2)); break; case 1: dofs.push_back(boundIt->rankObject.el->getDOF(0)); dofs.push_back(boundIt->rankObject.el->getDOF(2)); break; case 2: dofs.push_back(boundIt->rankObject.el->getDOF(0)); dofs.push_back(boundIt->rankObject.el->getDOF(1)); break; default: ERROR_EXIT("Should never happen!\n"); } #if 0 if (mpiRank == 3 && it->first == 2) { WorldVector x; mesh->getDofIndexCoords(dofs[0], feSpace, x); x.print(); mesh->getDofIndexCoords(dofs[1], feSpace, x); x.print(); } #endif for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt) if (find(dofsToSend.begin(), dofsToSend.end(), *dofIt) == dofsToSend.end()) dofsToSend.push_back(*dofIt); dofs.clear(); addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, dofs); addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, 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) { for (std::vector::iterator boundIt = it->second.begin(); boundIt != it->second.end(); ++boundIt) { #if 0 if (mpiRank == 2 && it->first == 3) std::cout << "GO ON OT BOUND: " << boundIt->rankObject.elIndex << "/" << boundIt->rankObject.ithObjAtBoundary << " with " << boundIt->neighbourObject.elIndex << "/" << boundIt->neighbourObject.ithObjAtBoundary << std::endl; #endif DofContainer dofs; DofContainer &dofsToRecv = recvDofs[it->first]; switch (boundIt->rankObject.ithObjAtBoundary) { case 0: if (boundIt->neighbourObject.ithObjAtBoundary == 0) { dofs.push_back(boundIt->rankObject.el->getDOF(2)); dofs.push_back(boundIt->rankObject.el->getDOF(1)); } else { dofs.push_back(boundIt->rankObject.el->getDOF(1)); dofs.push_back(boundIt->rankObject.el->getDOF(2)); } break; case 1: if (boundIt->neighbourObject.ithObjAtBoundary == 0) { dofs.push_back(boundIt->rankObject.el->getDOF(0)); dofs.push_back(boundIt->rankObject.el->getDOF(2)); } else { dofs.push_back(boundIt->rankObject.el->getDOF(2)); dofs.push_back(boundIt->rankObject.el->getDOF(0)); } break; case 2: dofs.push_back(boundIt->rankObject.el->getDOF(1)); dofs.push_back(boundIt->rankObject.el->getDOF(0)); break; default: ERROR_EXIT("Should never happen!\n"); } #if 0 if (mpiRank == 2 && it->first == 3) { WorldVector x; mesh->getDofIndexCoords(dofs[0], feSpace, x); x.print(); mesh->getDofIndexCoords(dofs[1], feSpace, x); x.print(); } #endif for (DofContainer::iterator dofIt = dofs.begin(); dofIt != dofs.end(); ++dofIt) { DofContainer::iterator eraseIt = find(rankDOFs.begin(), rankDOFs.end(), *dofIt); if (eraseIt != rankDOFs.end()) rankDOFs.erase(eraseIt); if (find(dofsToRecv.begin(), dofsToRecv.end(), *dofIt) == dofsToRecv.end()) dofsToRecv.push_back(*dofIt); } dofs.clear(); addAllEdgeDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, dofs); addAllVertexDOFs(boundIt->rankObject.el, boundIt->rankObject.ithObjAtBoundary, dofs); for (int i = static_cast(dofs.size()) - 1; i >= 0; i--) { TEST_EXIT_DBG(find(rankDOFs.begin(), rankDOFs.end(), dofs[i]) != rankDOFs.end()) ("Should never happen!\n"); 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. ==== int rstart = 0; 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::addAllVertexDOFs(Element *el, int ithEdge, DofContainer& dofs) { FUNCNAME("ParallelDomainBase::addAllVertexDOFs()"); switch (ithEdge) { case 0: if (el->getSecondChild() && el->getSecondChild()->getFirstChild()) { addAllVertexDOFs(el->getSecondChild()->getFirstChild(), 0, dofs); dofs.push_back(el->getSecondChild()->getFirstChild()->getDOF(2)); addAllVertexDOFs(el->getSecondChild()->getSecondChild(), 1, dofs); } break; case 1: if (el->getFirstChild() && el->getFirstChild()->getFirstChild()) { addAllVertexDOFs(el->getFirstChild()->getFirstChild(), 0, dofs); dofs.push_back(el->getFirstChild()->getFirstChild()->getDOF(2)); addAllVertexDOFs(el->getFirstChild()->getSecondChild(), 1, dofs); } break; case 2: if (el->getFirstChild()) { addAllVertexDOFs(el->getFirstChild(), 0, dofs); dofs.push_back(el->getFirstChild()->getDOF(2)); addAllVertexDOFs(el->getSecondChild(), 1, dofs); } break; default: ERROR_EXIT("Should never happen!\n"); } } void ParallelDomainBase::addAllEdgeDOFs(Element *el, int ithEdge, DofContainer& dofs) { FUNCNAME("ParallelDomainBase::addAllEdgeDOFs()"); bool addThisEdge = false; switch (ithEdge) { case 0: if (el->getSecondChild()) addAllEdgeDOFs(el->getSecondChild(), 2, dofs); else addThisEdge = true; break; case 1: if (el->getFirstChild()) addAllEdgeDOFs(el->getFirstChild(), 2, dofs); else addThisEdge = true; break; case 2: if (el->getFirstChild()) { addAllEdgeDOFs(el->getFirstChild(), 0, dofs); addAllEdgeDOFs(el->getSecondChild(), 1, dofs); } else { addThisEdge = true; } break; default: ERROR_EXIT("Should never happen!\n"); } if (addThisEdge) { ElementDofIterator elDofIter(feSpace, true); elDofIter.reset(el); do { if (elDofIter.getCurrentPos() == 1 && elDofIter.getCurrentElementPos() == ithEdge) dofs.push_back(elDofIter.getDofPtr()); } while(elDofIter.next()); } } void ParallelDomainBase::createDOFMemberInfo(DofToPartitions& partitionDofs, DofContainer& rankOwnedDofs, DofContainer& rankAllDofs, DofToRank& boundaryDofs) { partitionDofs.clear(); rankOwnedDofs.clear(); rankAllDofs.clear(); boundaryDofs.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()]); } 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) { // 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; } break; } } } sort(rankAllDofs.begin(), rankAllDofs.end(), cmpDofsByValue); sort(rankOwnedDofs.begin(), rankOwnedDofs.end(), cmpDofsByValue); } void ParallelDomainBase::DbgCreateElementMap(ElementIdxToDofs &elMap) { FUNCNAME("ParallelDomainBase::DbgCreateElementMap()"); TEST_EXIT(mesh->getDim() == 2)("Works only for 2d!\n"); elMap.clear(); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { Element *el = elInfo->getElement(); orderDOFs(el->getDOF(0), el->getDOF(1), el->getDOF(2), elMap[el->getIndex()]); elInfo = stack.traverseNext(elInfo); } } void ParallelDomainBase::DbgTestElementMap(ElementIdxToDofs &elMap) { FUNCNAME("ParallelDomainbase::DbgTestElementMap()"); TEST_EXIT(mesh->getDim() == 2)("Works only for 2d!\n"); DofContainer vec(3); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { Element *el = elInfo->getElement(); orderDOFs(el->getDOF(0), el->getDOF(1), el->getDOF(2), vec); for (int i = 0; i < 3; 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 < 3; 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 < 3; 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++) { #if 0 if (mpiRank == 3 && rankIt->first == 2) { std::cout << "MY BOUND: " << (rankIt->second)[i].rankObject.elIndex << "/" << (rankIt->second)[i].rankObject.ithObjAtBoundary << " with neig " << (rankIt->second)[i].neighbourObject.elIndex << "/" << (rankIt->second)[i].neighbourObject.ithObjAtBoundary << std::endl; } #endif buffer[i] = (rankIt->second)[i].rankObject.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++) { #if 0 if (mpiRank == 2 && rankIt->first == 3) { std::cout << "OT BOUND: " << (rankIt->second)[i].rankObject.elIndex << "/" << (rankIt->second)[i].rankObject.ithObjAtBoundary << " with neig " << (rankIt->second)[i].neighbourObject.elIndex << "/" << (rankIt->second)[i].neighbourObject.ithObjAtBoundary << std::endl; } #endif int elIndex1 = recvBuffers[bufCounter][i]; int elIndex2 = otherIntBoundary.boundary[rankIt->first][i].neighbourObject.elIndex; TEST_EXIT(elIndex1 == elIndex2)("Wrong element index at interior boundary!\n"); } delete [] recvBuffers[bufCounter++]; } } void ParallelDomainBase::DbgTestCommonDofs(bool printCoords) { FUNCNAME("ParallelDomainBase::DbgTestCommonDofs()"); #if 0 { TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_EVERY_EL_PREORDER); while (elInfo) { Element *el = elInfo->getElement(); for (int i = 0; i < 3; i++) { if (mpiRank == 2 && el->getDOF(i, 0) == 3) std::cout << "RANK 2 EL = " << el->getIndex() << " i = " << i << std::endl; if (mpiRank == 3 && el->getDOF(i, 0) == 4) std::cout << "RANK 3 EL = " << el->getIndex() << " i = " << i << std::endl; } elInfo = stack.traverseNext(elInfo); } } #endif // Maps to each neighbour rank an array of WorldVectors. This array contain 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); TEST_EXIT(b)("Cannot find DOF in mesh!\n"); #if 0 if (mpiRank == 3 && it->first == 2) { std::cout << "SEND COORDS: " << *dofIt << " = " << **dofIt << std::endl; coords.print(); } #endif 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); TEST_EXIT(b)("Cannot find DOF in mesh!\n"); recvCoords[it->first].push_back(coords); #if 0 if (mpiRank == 2 && it->first == 3) { std::cout << "RECV COORDS: " << *dofIt << " = " << **dofIt << std::endl; coords.print(); } #endif } } 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] << " , "; std::cout << ") received coords ("; for (int k = 0; k < dimOfWorld; k++) std::cout << recvCoordsBuffer[it->first][i * dimOfWorld + k] << " , "; std::cout << ")" << std::endl; } TEST_EXIT(c)("Wrong DOFs in rank %d!\n", mpiRank); } } delete [] recvCoordsBuffer[it->first]; } } Flag ParallelDomainBase::oneIteration(AdaptInfo *adaptInfo, Flag toDo) { FUNCNAME("ParallelDomainBase::oneIteration()"); Flag flag = dynamic_cast(iterationIF)-> 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; } }