diff --git a/AMDiS/src/InteriorBoundary.cc b/AMDiS/src/InteriorBoundary.cc index 869f3705f87b3c4a1c870b5bdd2ee13ae86adfe0..830d1efcdad159f69740c1362cbbc60d79ccf3a1 100644 --- a/AMDiS/src/InteriorBoundary.cc +++ b/AMDiS/src/InteriorBoundary.cc @@ -18,7 +18,9 @@ namespace AMDiS { switch (feSpace->getMesh()->getDim()) { case 2: - ERROR_EXIT("Not yet implemented!\n"); + if (ithObj == 2 || ithObj != otherBound.ithObj) + otherMode = true; + break; case 3: if (ithObj == 2 || ithObj == 3) { diff --git a/AMDiS/src/InteriorBoundary.h b/AMDiS/src/InteriorBoundary.h index ae85f77b869304f8f962973a258e37fe4b392a2d..3dc8d9572a751a65cb2867e73a95583e67038ec0 100644 --- a/AMDiS/src/InteriorBoundary.h +++ b/AMDiS/src/InteriorBoundary.h @@ -116,6 +116,82 @@ namespace AMDiS { public: typedef std::map<int, std::vector<AtomicBoundary> > RankToBoundMap; + /// Iterator for the interior boundary object. + class iterator { + public: + iterator(InteriorBoundary &b) + : bound(b) + { + reset(); + } + + /// Set the iterator to the first position. + void reset() + { + mapIt = bound.boundary.begin(); + nextNonempty(); + + if (mapIt != bound.boundary.end()) + vecIt = mapIt->second.begin(); + } + + /// Test if iterator is at the final position. + bool end() const + { + return (mapIt == bound.boundary.end()); + } + + /// Move iterator to the next position. + void operator++() + { + ++vecIt; + if (vecIt == mapIt->second.end()) { + ++mapIt; + nextNonempty(); + + if (mapIt != bound.boundary.end()) + vecIt = mapIt->second.begin(); + } + } + + void nextRank() + { + ++mapIt; + nextNonempty(); + + if (mapIt != bound.boundary.end()) + vecIt = mapIt->second.begin(); + } + + int getRank() + { + return mapIt->first; + } + + AtomicBoundary& getBound() + { + return *vecIt; + } + + protected: + + inline void nextNonempty() + { + while (mapIt->second.size() == 0) { + ++mapIt; + if (mapIt == bound.boundary.end()) + return; + } + } + + protected: + RankToBoundMap::iterator mapIt; + + std::vector<AtomicBoundary>::iterator vecIt; + + InteriorBoundary &bound; + }; + public: InteriorBoundary() {} diff --git a/AMDiS/src/ParallelDomainBase.cc b/AMDiS/src/ParallelDomainBase.cc index 0cca1611fb7e80b8bfa8751800dff625a2f15ea8..f4327b943f3e5e4d26cef12ad8995d10504a4507 100644 --- a/AMDiS/src/ParallelDomainBase.cc +++ b/AMDiS/src/ParallelDomainBase.cc @@ -730,6 +730,8 @@ namespace AMDiS { { FUNCNAME("ParallelDomainBase::checkMeshChange()"); + int dim = mesh->getDim(); + // === If mesh has not been changed on all ranks, return. === int recvAllValues = 0; @@ -748,19 +750,15 @@ namespace AMDiS { // To check the interior boundaries, the ownership of the boundaries is not // important. Therefore, we add all boundaries to one boundary container. RankToBoundMap allBound; - for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin(); - it != myIntBoundary.boundary.end(); ++it) - for (std::vector<AtomicBoundary>::iterator bit = it->second.begin(); - bit != it->second.end(); ++bit) - if (bit->rankObj.subObj == FACE) - allBound[it->first].push_back(*bit); - - for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin(); - it != otherIntBoundary.boundary.end(); ++it) - for (std::vector<AtomicBoundary>::iterator bit = it->second.begin(); - bit != it->second.end(); ++bit) - if (bit->rankObj.subObj == FACE) - allBound[it->first].push_back(*bit); + + for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) + if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim)) + allBound[it.getRank()].push_back(it.getBound()); + + for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) + if (it.getBound().rankObj.subObj == INDEX_OF_DIM(dim - 1, dim)) + allBound[it.getRank()].push_back(it.getBound()); + // === Check the boundaries and adapt mesh if necessary. === @@ -771,7 +769,6 @@ namespace AMDiS { int sendValue = static_cast<int>(!meshChanged); recvAllValues = 0; mpiComm.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM); - MSG("LOOP\n"); } while (recvAllValues != 0); @@ -805,8 +802,6 @@ namespace AMDiS { { FUNCNAME("ParallelDomainBase::checkAndAdaptBoundary()"); - MSG("CHECK_AND_ADAPT 1!\n"); - // === Create mesh structure codes for all ranks boundary elements. === std::map<int, MeshCodeVec> sendCodes; @@ -820,21 +815,14 @@ namespace AMDiS { } } - MSG("CHECK_AND_ADAPT 2!\n"); - StdMpi<MeshCodeVec> stdMpi(mpiComm, true); stdMpi.send(sendCodes); stdMpi.recv(allBound); stdMpi.startCommunication<unsigned long int>(MPI_UNSIGNED_LONG); - MSG("CHECK_AND_ADAPT 3!\n"); - - clock_t first = clock(); - // === Compare received mesh structure codes. === bool meshFitTogether = true; - int superCounter = 0; for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) { @@ -849,8 +837,6 @@ namespace AMDiS { if (elCode.getCode() != recvCodes[i].getCode()) { TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n"); - superCounter++; - clock_t fff = clock(); int refCounter = 0; bool b = fitElementToMeshCode(recvCodes[i], @@ -858,12 +844,6 @@ namespace AMDiS { boundIt->rankObj.ithObj, boundIt->rankObj.elType, refCounter); - MSG("SIZE OF ELINFO3d = %d\n", sizeof(ElInfo3d)); - MSG("Code length %d with %d refs needed %.5f seconds\n", - recvCodes[i].getNumElements(), - refCounter, - TIME_USED(fff, clock())); - if (b) meshFitTogether = false; } @@ -872,9 +852,6 @@ namespace AMDiS { } } - MSG("time for %d needed %.5f seconds\n", superCounter, TIME_USED(first, clock())); - MSG("CHECK_AND_ADAPT 4!\n"); - return meshFitTogether; } @@ -1370,77 +1347,69 @@ namespace AMDiS { std::set<DegreeOfFreedom> rankBoundaryDofs; // First, traverse the rank owned elements af the interior boundaries. - for (RankToBoundMap::iterator rankIt = myIntBoundary.boundary.begin(); - rankIt != myIntBoundary.boundary.end(); ++rankIt) { - for (std::vector<AtomicBoundary>::iterator boundIt = rankIt->second.begin(); - boundIt != rankIt->second.end(); ++boundIt) { - Element *el = boundIt->rankObj.el; - basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); - - if (dim == 3) { - for (int j = 0; j < 3; j++) { - int edgeNo = el->getEdgeOfFace(boundIt->rankObj.ithObj, j); - DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)]; - DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)]; - GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); - - // If the edge at rank's interior boundarie has a higher owner rank, than - // we have to remove this edge from the corresponding boundary element. - // Otherwise, it is part of the interior boundary and we add it to the set - // rankBoundaryEdges. - if (edgeOwner[edge] > mpiRank) - boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo)); - else - rankBoundaryEdges.insert(edge); - } - } + for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) { + Element *el = it.getBound().rankObj.el; + basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); + if (dim == 3) { for (int j = 0; j < 3; j++) { - int dofNo = el->getVertexOfPosition(FACE, boundIt->rankObj.ithObj, j); - DegreeOfFreedom dof = localIndices[dofNo]; - - if (dofOwner[dof] > mpiRank) - boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo)); + int edgeNo = el->getEdgeOfFace(it.getBound().rankObj.ithObj, j); + DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)]; + DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)]; + GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); + + // If the edge at rank's interior boundarie has a higher owner rank, than + // we have to remove this edge from the corresponding boundary element. + // Otherwise, it is part of the interior boundary and we add it to the set + // rankBoundaryEdges. + if (edgeOwner[edge] > mpiRank) + it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo)); else - rankBoundaryDofs.insert(dof); + rankBoundaryEdges.insert(edge); } } + + for (int j = 0; j < 3; j++) { + int dofNo = el->getVertexOfPosition(FACE, it.getBound().rankObj.ithObj, j); + DegreeOfFreedom dof = localIndices[dofNo]; + + if (dofOwner[dof] > mpiRank) + it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo)); + else + rankBoundaryDofs.insert(dof); + } } - // Now, do the same with all other elements at the interio boundaries. - for (RankToBoundMap::iterator rankIt = otherIntBoundary.boundary.begin(); - rankIt != otherIntBoundary.boundary.end(); ++rankIt) { - for (std::vector<AtomicBoundary>::iterator boundIt = rankIt->second.begin(); - boundIt != rankIt->second.end(); ++boundIt) { - Element *el = boundIt->rankObj.el; - basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); - + // Now, do the same with all other elements at the interior boundaries. + for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) { + Element *el = it.getBound().rankObj.el; + basFcts->getLocalIndices(el, feSpace->getAdmin(), localIndices); + + if (dim == 3) { for (int j = 0; j < 3; j++) { - int edgeNo = el->getEdgeOfFace(boundIt->rankObj.ithObj, j); + int edgeNo = el->getEdgeOfFace(it.getBound().rankObj.ithObj, j); DegreeOfFreedom dof0 = localIndices[el->getVertexOfEdge(edgeNo, 0)]; DegreeOfFreedom dof1 = localIndices[el->getVertexOfEdge(edgeNo, 1)]; GlobalEdge edge = std::make_pair(min(dof0, dof1), max(dof0, dof1)); - - if (edgeOwner[edge] > rankIt->first) - boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo)); - else - rankBoundaryEdges.insert(edge); - - } - - for (int j = 0; j < 3; j++) { - int dofNo = el->getVertexOfPosition(FACE, boundIt->rankObj.ithObj, j); - DegreeOfFreedom dof = localIndices[dofNo]; - if (dofOwner[dof] > rankIt->first) - boundIt->rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo)); + if (edgeOwner[edge] > it.getRank()) + it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(EDGE, edgeNo)); else - rankBoundaryDofs.insert(dof); - } + rankBoundaryEdges.insert(edge); + } } + + for (int j = 0; j < 3; j++) { + int dofNo = el->getVertexOfPosition(FACE, it.getBound().rankObj.ithObj, j); + DegreeOfFreedom dof = localIndices[dofNo]; + + if (dofOwner[dof] > it.getRank()) + it.getBound().rankObj.excludedSubstructures.push_back(std::make_pair(VERTEX, dofNo)); + else + rankBoundaryDofs.insert(dof); + } } - // === Create the new interior boundaries consisting only of edges. This === // === boundaries are created on that ranks, which do not own the boundary === // === but are on the other side of the edge. Than, theses ranks inform the === @@ -1894,45 +1863,29 @@ namespace AMDiS { } } - - for (RankToBoundMap::iterator it = myIntBoundary.boundary.begin(); - it != myIntBoundary.boundary.end(); ++it) { - - DofContainer &dofsToSend = sendDofs[it->first]; - - for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin(); - boundIt != it->second.end(); ++boundIt) { - DofContainer dofs; - boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, dofs); - boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj, dofs); - - for (int i = 0; i < static_cast<int>(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<AtomicBoundary>::iterator boundIt = it->second.begin(); - boundIt != it->second.end(); ++boundIt) { - DofContainer dofs; - boundIt->rankObj.el->getNonVertexDofs(feSpace, boundIt->rankObj, dofs); - boundIt->rankObj.el->getVertexDofs(feSpace, boundIt->rankObj, dofs); - - for (int i = 0; i < static_cast<int>(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]); - } - } + for (InteriorBoundary::iterator it(myIntBoundary); !it.end(); ++it) { + DofContainer dofs; + it.getBound().rankObj.el->getVertexDofs(feSpace, it.getBound().rankObj, dofs); + it.getBound().rankObj.el->getNonVertexDofs(feSpace, it.getBound().rankObj, dofs); + for (int i = 0; i < static_cast<int>(dofs.size()); i++) + sendDofs[it.getRank()].push_back(dofs[i]); } + for (InteriorBoundary::iterator it(otherIntBoundary); !it.end(); ++it) { + DofContainer dofs; + it.getBound().rankObj.el->getNonVertexDofs(feSpace, it.getBound().rankObj, dofs); + it.getBound().rankObj.el->getVertexDofs(feSpace, it.getBound().rankObj, dofs); + + for (int i = 0; i < static_cast<int>(dofs.size()); i++) { + DofContainer::iterator eraseIt = + find(rankDofs.begin(), rankDofs.end(), dofs[i]); + if (eraseIt != rankDofs.end()) + rankDofs.erase(eraseIt); + + recvDofs[it.getRank()].push_back(dofs[i]); + } + } + nRankDofs = rankDofs.size(); // === Get starting position for global rank dof ordering. ====