diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc index 0e715f6ea0764e8ab6df87186bade2cdcf9c4178..fb3ac5b4b328b5283dc7403b2184c6e611d022e1 100644 --- a/AMDiS/src/Debug.cc +++ b/AMDiS/src/Debug.cc @@ -768,6 +768,30 @@ namespace AMDiS { } + void testDofsByCoords(FiniteElemSpace *feSpace, + DofContainer &dofs0, DofContainer &dofs1) + { + FUNCNAME("debug::testDofsByCoords()"); + + TEST_EXIT(dofs0.size() == dofs1.size()) + ("The dof containers have different sizes %d %d!\n", + dofs0.size(), dofs1.size()); + + DOFVector<WorldVector<double> > coords(feSpace, "dofCorrds"); + feSpace->getMesh()->getDofIndexCoords(feSpace, coords); + + for (unsigned int i = 0; i < dofs0.size(); i++) { + WorldVector<double> tmp = coords[*(dofs0[i])]; + tmp -= coords[*(dofs1[i])]; + + TEST_EXIT(norm(tmp) < 1e-13) + ("DOFs %d and %d (i = %d) have different coords!\n", + *(dofs0[i]), *(dofs1[i]), i); + } + + } + + } // namespace debug } // namespace AMDiS diff --git a/AMDiS/src/Debug.h b/AMDiS/src/Debug.h index 458027214c813375dda671ca660e0f312c27b850..9d0d141b8b2a3f4a2741982726b7684dd7311ccf 100644 --- a/AMDiS/src/Debug.h +++ b/AMDiS/src/Debug.h @@ -193,6 +193,9 @@ namespace AMDiS { const DegreeOfFreedom* dof3, DofContainer &vec); + void testDofsByCoords(FiniteElemSpace *feSpace, + DofContainer &dofs0, DofContainer &dofs1); + } } diff --git a/AMDiS/src/parallel/ElementObjectData.cc b/AMDiS/src/parallel/ElementObjectData.cc index 065f1224ca8ddd9a0de39fd59e5f17c8622a9b4a..18c072be4de528ba44ab5fc1e57cdc51c00cfc68 100644 --- a/AMDiS/src/parallel/ElementObjectData.cc +++ b/AMDiS/src/parallel/ElementObjectData.cc @@ -97,6 +97,32 @@ namespace AMDiS { + nSize = vertexLocalMap.size(); + SerUtil::serialize(out, nSize); + for (std::map<ElementObjectData, DegreeOfFreedom>::iterator it = vertexLocalMap.begin(); + it != vertexLocalMap.end(); ++it) { + it->first.serialize(out); + SerUtil::serialize(out, it->second); + } + + nSize = edgeLocalMap.size(); + SerUtil::serialize(out, nSize); + for (std::map<ElementObjectData, DofEdge>::iterator it = edgeLocalMap.begin(); + it != edgeLocalMap.end(); ++it) { + it->first.serialize(out); + SerUtil::serialize(out, it->second); + } + + nSize = faceLocalMap.size(); + SerUtil::serialize(out, nSize); + for (std::map<ElementObjectData, DofFace>::iterator it = faceLocalMap.begin(); + it != faceLocalMap.end(); ++it) { + it->first.serialize(out); + SerUtil::serialize(out, it->second); + } + + + SerUtil::serialize(out, vertexOwner); SerUtil::serialize(out, edgeOwner); SerUtil::serialize(out, faceOwner); @@ -165,6 +191,38 @@ namespace AMDiS { } + + SerUtil::deserialize(in, nSize); + vertexLocalMap.clear(); + for (int i = 0; i < nSize; i++) { + ElementObjectData data; + DegreeOfFreedom dof; + data.deserialize(in); + SerUtil::deserialize(in, dof); + vertexLocalMap[data] = dof; + } + + SerUtil::deserialize(in, nSize); + edgeLocalMap.clear(); + for (int i = 0; i < nSize; i++) { + ElementObjectData data; + DofEdge edge; + data.deserialize(in); + SerUtil::deserialize(in, edge); + edgeLocalMap[data] = edge; + } + + SerUtil::deserialize(in, nSize); + faceLocalMap.clear(); + for (int i = 0; i < nSize; i++) { + ElementObjectData data; + DofFace face; + data.deserialize(in); + SerUtil::deserialize(in, face); + faceLocalMap[data] = face; + } + + SerUtil::deserialize(in, vertexOwner); SerUtil::deserialize(in, edgeOwner); SerUtil::deserialize(in, faceOwner); diff --git a/AMDiS/src/parallel/ElementObjectData.h b/AMDiS/src/parallel/ElementObjectData.h index 7975644a3d79c0c7167c437a0ba9764b4733cde6..883c9a455675c90b3a801cb10c32c7f1450d916e 100644 --- a/AMDiS/src/parallel/ElementObjectData.h +++ b/AMDiS/src/parallel/ElementObjectData.h @@ -48,7 +48,7 @@ namespace AMDiS { BoundaryType boundaryType; - void serialize(std::ostream &out) + void serialize(std::ostream &out) const { SerUtil::serialize(out, elIndex); SerUtil::serialize(out, ithObject); @@ -63,6 +63,20 @@ namespace AMDiS { SerUtil::deserialize(in, boundaryType); } + + bool operator==(ElementObjectData& cmp) const + { + return (elIndex == cmp.elIndex && + ithObject == cmp.ithObject && + boundaryType == cmp.boundaryType); + } + + bool operator<(const ElementObjectData& rhs) const + { + return (elIndex < rhs.elIndex || + (elIndex == rhs.elIndex && + ithObject < rhs.ithObject)); + } }; @@ -73,26 +87,56 @@ namespace AMDiS { : iterGeoPos(CENTER) {} - void addVertex(DegreeOfFreedom vertex, - int elIndex, int ith, BoundaryType bound = INTERIOR) + + void addVertex(Element *el, int ith, BoundaryType bound = INTERIOR) { - vertexElements[vertex].push_back(ElementObjectData(elIndex, ith, bound)); + DegreeOfFreedom vertex = el->getDof(ith, 0); + int elIndex = el->getIndex(); + ElementObjectData elObj(elIndex, ith, bound); + + vertexElements[vertex].push_back(elObj); + vertexLocalMap[elObj] = vertex; } - void addEdge(DofEdge edge, - int elIndex, int ith, BoundaryType bound = INTERIOR) + + void addEdge(Element *el, int ith, BoundaryType bound = INTERIOR) { - edgeElements[edge].push_back(ElementObjectData(elIndex, ith, bound)); + DofEdge edge = el->getEdge(ith); + int elIndex = el->getIndex(); + ElementObjectData elObj(elIndex, ith, bound); + + edgeElements[edge].push_back(elObj); + edgeLocalMap[elObj] = edge; + } + + + void addFace(Element *el, int ith, BoundaryType bound = INTERIOR) + { + DofFace face = el->getFace(ith); + int elIndex = el->getIndex(); + ElementObjectData elObj(elIndex, ith, bound); + + faceElements[face].push_back(elObj); + faceLocalMap[elObj] = face; } - void addFace(DofFace face, - int elIndex, int ith, BoundaryType bound = INTERIOR) + + void addElement(Element *el, BoundaryType bound = INTERIOR) { - faceElements[face].push_back(ElementObjectData(elIndex, ith, bound)); + for (int i = 0; i < el->getGeo(VERTEX); i++) + addVertex(el, i); + + for (int i = 0; i < el->getGeo(EDGE); i++) + addEdge(el, i); + + for (int i = 0; i < el->getGeo(FACE); i++) + addFace(el, i); } + void createRankData(std::map<int, int>& macroElementRankMap); + bool iterate(GeoIndex pos) { if (iterGeoPos == CENTER) { @@ -258,6 +302,21 @@ namespace AMDiS { return faceInRank[face]; } + DegreeOfFreedom getVertexLocalMap(ElementObjectData &data) + { + return vertexLocalMap[data]; + } + + DofEdge getEdgeLocalMap(ElementObjectData &data) + { + return edgeLocalMap[data]; + } + + DofFace getFaceLocalMap(ElementObjectData &data) + { + return faceLocalMap[data]; + } + void serialize(std::ostream &out); void deserialize(std::istream &in); @@ -276,14 +335,22 @@ namespace AMDiS { std::map<DofEdge, std::vector<ElementObjectData> > edgeElements; std::map<DofFace, std::vector<ElementObjectData> > faceElements; + + std::map<ElementObjectData, DegreeOfFreedom> vertexLocalMap; + std::map<ElementObjectData, DofEdge> edgeLocalMap; + std::map<ElementObjectData, DofFace> faceLocalMap; + + std::map<DegreeOfFreedom, int> vertexOwner; std::map<DofEdge, int> edgeOwner; std::map<DofFace, int> faceOwner; + std::map<DegreeOfFreedom, std::map<int, ElementObjectData> > vertexInRank; std::map<DofEdge, std::map<int, ElementObjectData> > edgeInRank; std::map<DofFace, std::map<int, ElementObjectData> > faceInRank; + std::map<DegreeOfFreedom, std::map<int, ElementObjectData> >::iterator vertexIter; std::map<DofEdge, std::map<int, ElementObjectData> >::iterator edgeIter; std::map<DofFace, std::map<int, ElementObjectData> >::iterator faceIter; diff --git a/AMDiS/src/parallel/InteriorBoundary.h b/AMDiS/src/parallel/InteriorBoundary.h index ccbdbc5b5e690bca8e5919264bbf225bc43b248e..ea12b82f7fb8322cfb3d43e6f964ea3fa9f67a9b 100644 --- a/AMDiS/src/parallel/InteriorBoundary.h +++ b/AMDiS/src/parallel/InteriorBoundary.h @@ -44,13 +44,13 @@ namespace AMDiS { excludedSubstructures(0) {} - BoundaryObject(Element *e, int eType, GeoIndex sObj, int iObj) + BoundaryObject(Element *e, int eType, GeoIndex sObj, int iObj, bool rMode = false) : el(e), elIndex(e->getIndex()), elType(eType), subObj(sObj), ithObj(iObj), - reverseMode(false), + reverseMode(rMode), excludedSubstructures(0) {} diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index 4f21cd562bf0c2be00519f502a6c31afa07297c4..a169e2b6b5c4dff89dd90631da576700a221a4e0 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -1,4 +1,4 @@ - +// // Software License for AMDiS // // Copyright (c) 2010 Dresden University of Technology @@ -68,6 +68,7 @@ namespace AMDiS { deserialized(false), writeSerializationFile(false), repartitioningAllowed(false), + repartitionIthChange(20), nTimestepsAfterLastRepartitioning(0), repartCounter(0), debugOutputDir(""), @@ -83,7 +84,9 @@ namespace AMDiS { GET_PARAMETER(0, name + "->repartitioning", "%d", &tmp); repartitioningAllowed = (tmp > 0); - GET_PARAMETER(0, name + "->debug output dir", &debugOutputDir); + GET_PARAMETER(0, name + "->debug output dir", &debugOutputDir); + + GET_PARAMETER(0, name + "->repartition ith change", "%d", &repartitionIthChange); } @@ -581,7 +584,7 @@ namespace AMDiS { nTimestepsAfterLastRepartitioning++; if (repartitioningAllowed) { - if (nTimestepsAfterLastRepartitioning >= 20) { + if (nTimestepsAfterLastRepartitioning >= repartitionIthChange) { repartitionMesh(); nTimestepsAfterLastRepartitioning = 0; } @@ -1001,7 +1004,9 @@ namespace AMDiS { return; #if (DEBUG != 0) - if (repartCounter == 0) { + ParallelDebug::testDoubleDofs(mesh); + + if (repartCounter == 0) { std::stringstream oss; oss << debugOutputDir << "partitioning-" << repartCounter << ".vtu"; DOFVector<double> tmpa(feSpace, "tmp"); @@ -1191,8 +1196,8 @@ namespace AMDiS { // === Remove double DOFs. === - MeshManipulation meshManipulation(mesh); - meshManipulation.deleteDoubleDofs(newMacroEl); + MeshManipulation meshManipulation(feSpace); + meshManipulation.deleteDoubleDofs(newMacroEl, elObjects); mesh->dofCompress(); @@ -1207,21 +1212,24 @@ namespace AMDiS { MSG("USED-SIZE B: %d\n", mesh->getDofAdmin(0).getUsedDofs()); ParallelDebug::testAllElements(*this); + ParallelDebug::testDoubleDofs(mesh); #endif partitioner->fillCoarsePartitionVec(&partitionVec); updateInteriorBoundaryInfo(); +#if (DEBUG != 0) + ParallelDebug::printBoundaryInfo(*this); +#endif + updateLocalGlobalNumbering(); - #if (DEBUG != 0) +#if (DEBUG != 0) MSG("AMDiS runs in debug mode, so make some test ...\n"); ParallelDebug::testAllElements(*this); ParallelDebug::testInteriorBoundary(*this); - ParallelDebug::testCommonDofs(*this, true); - ParallelDebug::testGlobalIndexByCoords(*this); debug::writeMesh(feSpace, -1, debugOutputDir + "macro_mesh"); @@ -1273,14 +1281,7 @@ namespace AMDiS { // === Add all sub object of the element to the variable elObjects. === - for (int i = 0; i < el->getGeo(VERTEX); i++) - elObjects.addVertex(el->getDof(i, 0), el->getIndex(), i); - - for (int i = 0; i < el->getGeo(EDGE); i++) - elObjects.addEdge(el->getEdge(i), el->getIndex(), i); - - for (int i = 0; i < el->getGeo(FACE); i++) - elObjects.addFace(el->getFace(i), el->getIndex(), i); + elObjects.addElement(el); // === Get periodic boundary information. === @@ -1801,7 +1802,8 @@ namespace AMDiS { // === Update dof admins due to new number of dofs. === lastMeshChangeIndex = mesh->getChangeIndex(); - + + #if (DEBUG != 0) std::stringstream oss; oss << debugOutputDir << "elementIndex-" << mpiRank << ".vtu"; diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h index 3eac43f564f89021f743e9f608f3b292497a6975..d4597e8ce3405c653508b209278234062bd24201 100644 --- a/AMDiS/src/parallel/MeshDistributor.h +++ b/AMDiS/src/parallel/MeshDistributor.h @@ -568,6 +568,8 @@ namespace AMDiS { /// If true, it is possible to repartition the mesh during computations. bool repartitioningAllowed; + int repartitionIthChange; + int nTimestepsAfterLastRepartitioning; int repartCounter; diff --git a/AMDiS/src/parallel/MeshManipulation.cc b/AMDiS/src/parallel/MeshManipulation.cc index f894f5515b3425e074f23ae04a85e005ef36ff1b..ebf60fa7f7d920d884b2b9129fd5ac55243f1808 100644 --- a/AMDiS/src/parallel/MeshManipulation.cc +++ b/AMDiS/src/parallel/MeshManipulation.cc @@ -12,118 +12,118 @@ #include "parallel/MeshManipulation.h" #include "Mesh.h" +#include "BasisFunction.h" #include "Traverse.h" #include "Debug.h" namespace AMDiS { - void MeshManipulation::deleteDoubleDofs(std::set<MacroElement*>& newMacroEl) + void MeshManipulation::deleteDoubleDofs(std::set<MacroElement*>& newMacroEl, + ElementObjects &objects) { - std::map<int, MacroElement*> leafInMacroEl; - + FUNCNAME("MeshManipulation::deleteDoubleDofs()"); + + TEST_EXIT(mesh->getDim() == 2)("Not yet supported for dim != 2!\n"); + + std::map<int, MacroElement*> macroIndexMap; + for (std::set<MacroElement*>::iterator it = newMacroEl.begin(); + it != newMacroEl.end(); ++it) + macroIndexMap[(*it)->getIndex()] = *it; + + std::set<int> macrosProcessed; + std::map<const DegreeOfFreedom*, const DegreeOfFreedom*> mapDelDofs; + + TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); + ElInfo *elInfo = stack.traverseFirst(mesh, 0, Mesh::CALL_EL_LEVEL); while (elInfo) { - leafInMacroEl[elInfo->getElement()->getIndex()] = elInfo->getMacroElement(); + if (newMacroEl.count(elInfo->getMacroElement()) == 0) { + int index = elInfo->getMacroElement()->getIndex(); + + macrosProcessed.insert(index); + macroIndexMap[index] = elInfo->getMacroElement(); + } + elInfo = stack.traverseNext(elInfo); } - deleteDoubleDofs(leafInMacroEl, newMacroEl, 0); - deleteDoubleDofs(leafInMacroEl, newMacroEl, 1); - } + for (std::set<MacroElement*>::iterator it = newMacroEl.begin(); + it != newMacroEl.end(); ++it) { - void MeshManipulation::deleteDoubleDofs(std::map<int, MacroElement*>& leafInMacroEl, - std::set<MacroElement*>& newMacroEl, - int mode) - { - FUNCNAME("MeshManipulation::deleteDoubleDofs()"); + for (int i = 0; i < mesh->getGeo(VERTEX); i++) { + ElementObjectData elObj((*it)->getIndex(), i); + DegreeOfFreedom vertex = objects.getVertexLocalMap(elObj); + std::vector<ElementObjectData> &vertexEl = objects.getElements(vertex); - std::map<const DegreeOfFreedom*, const DegreeOfFreedom*> mapDelDofs; - std::map<const DegreeOfFreedom*, int> mapDofsMacro; + for (std::vector<ElementObjectData>::iterator elIt = vertexEl.begin(); + elIt != vertexEl.end(); ++elIt) { + if (elIt->elIndex == (*it)->getIndex()) + continue; - TraverseStack stack; - ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_NEIGH); - while (elInfo) { - Element *el = elInfo->getElement(); - - for (int i = 0; i < mesh->getGeo(NEIGH); i++) { - - Element *neigh = elInfo->getNeighbour(i); - - if (!neigh) - continue; - - TEST_EXIT_DBG(leafInMacroEl.count(el->getIndex()) == 1) - ("Should not happen!\n"); - TEST_EXIT_DBG(leafInMacroEl.count(neigh->getIndex()) == 1) - ("No macro element for el %d, that is %d-ith neigbour of element %d!\n", - neigh->getIndex(), i, el->getIndex()); - - int elMacroIndex = leafInMacroEl[el->getIndex()]->getIndex(); - int neighMacroIndex = leafInMacroEl[neigh->getIndex()]->getIndex(); - - if (elMacroIndex == neighMacroIndex) - continue; - - if ((mode == 0 && - newMacroEl.count(leafInMacroEl[el->getIndex()]) == 1 && - newMacroEl.count(leafInMacroEl[neigh->getIndex()]) == 1 && - elMacroIndex > neighMacroIndex) || - (mode == 1 && - newMacroEl.count(leafInMacroEl[el->getIndex()]) == 0 && - newMacroEl.count(leafInMacroEl[neigh->getIndex()]) == 1)) { - - if (el->getEdge(i) != neigh->getEdge(elInfo->getOppVertex(i))) { - std::vector<int> elEdgeIndex(2); - elEdgeIndex[0] = el->getVertexOfEdge(i, 0); - elEdgeIndex[1] = el->getVertexOfEdge(i, 1); - - std::vector<int> neighEdgeIndex(2); - neighEdgeIndex[0] = neigh->getVertexOfEdge(elInfo->getOppVertex(i), 0); - neighEdgeIndex[1] = neigh->getVertexOfEdge(elInfo->getOppVertex(i), 1); - - std::vector<DegreeOfFreedom> elEdgeDof(2); - elEdgeDof[0] = el->getDof(elEdgeIndex[0], 0); - elEdgeDof[1] = el->getDof(elEdgeIndex[1], 0); - - std::vector<DegreeOfFreedom> neighEdgeDof(2); - neighEdgeDof[0] = neigh->getDof(neighEdgeIndex[0], 0); - neighEdgeDof[1] = neigh->getDof(neighEdgeIndex[1], 0); - - if ((elEdgeDof[0] < elEdgeDof[1] && neighEdgeDof[0] > neighEdgeDof[1]) || - (elEdgeDof[0] > elEdgeDof[1] && neighEdgeDof[0] < neighEdgeDof[1])) { - std::swap(neighEdgeIndex[0], neighEdgeIndex[1]); - std::swap(neighEdgeDof[0], neighEdgeDof[1]); - } - - - for (int j = 0; j < 2; j++) { - - if (neighEdgeDof[j] != elEdgeDof[j]) { - const DegreeOfFreedom *delDof = neigh->getDof(neighEdgeIndex[j]); - const DegreeOfFreedom *replaceDof = el->getDof(elEdgeIndex[j]); - - if (mapDelDofs.count(delDof) == 1 && mapDelDofs[delDof] != replaceDof) { - if (mapDofsMacro[mapDelDofs[delDof]] > elMacroIndex) { - mapDelDofs[replaceDof] = mapDelDofs[delDof]; - } else { - mapDelDofs[mapDelDofs[delDof]] = replaceDof; - mapDelDofs[delDof] = replaceDof; - mapDofsMacro[replaceDof] = elMacroIndex; - } - } else { - mapDelDofs[delDof] = replaceDof; - mapDofsMacro[replaceDof] = elMacroIndex; - } - } - } + if (macrosProcessed.count(elIt->elIndex) == 1) { + TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex) == 1) + ("Should not happen!\n"); + + Element *el0 = (*it)->getElement(); + Element *el1 = macroIndexMap[elIt->elIndex]->getElement(); + + const DegreeOfFreedom *dof0 = el0->getDof(i); + const DegreeOfFreedom *dof1 = el1->getDof(elIt->ithObject); + + mapDelDofs[dof0] = dof1; + + break; + } + } + } + + for (int i = 0; i < mesh->getGeo(EDGE); i++) { + ElementObjectData elObj((*it)->getIndex(), i); + DofEdge edge = objects.getEdgeLocalMap(elObj); + std::vector<ElementObjectData> &edgeEl = objects.getElements(edge); + + for (std::vector<ElementObjectData>::iterator elIt = edgeEl.begin(); + elIt != edgeEl.end(); ++elIt) { + if (elIt->elIndex == (*it)->getIndex()) + continue; + + if (macrosProcessed.count(elIt->elIndex) == 1) { + TEST_EXIT_DBG(macroIndexMap.count(elIt->elIndex) == 1) + ("Should not happen!\n"); + + TEST_EXIT(mesh->getDim() == 2)("In 3D set correct reverse mode!\n"); + + TEST_EXIT(feSpace->getBasisFcts()->getDegree() == 1) + ("Not yet implemented!\n"); + + Element *el0 = (*it)->getElement(); + Element *el1 = macroIndexMap[elIt->elIndex]->getElement(); + + BoundaryObject b0(el0, 0, EDGE, i, true); + BoundaryObject b1(el1, 0, EDGE, elIt->ithObject, false); + + DofContainer dofs0, dofs1; + + el0->getVertexDofs(feSpace, b0, dofs0); + el1->getVertexDofs(feSpace, b1, dofs1); + +#if (DEBUG != 0) + debug::testDofsByCoords(feSpace, dofs0, dofs1); +#endif + + for (unsigned int i = 0; i < dofs0.size(); i++) + mapDelDofs[dofs0[i]] = dofs1[i]; + + break; } } } + TEST_EXIT(mesh->getDim() == 2)("Add face traverse here!\n"); - elInfo = stack.traverseNext(elInfo); + + macrosProcessed.insert((*it)->getIndex()); } @@ -134,6 +134,9 @@ namespace AMDiS { it != mapDelDofs.end(); ++it) { std::map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator findIt = mapDelDofs.find(it->second); if (findIt != mapDelDofs.end()) { + TEST_EXIT_DBG(it->first != findIt->second) + ("Cycle %d -> %d and %d -> %d! Should not happen!\n", + *(it->first), *(it->second), *(findIt->first), *(findIt->second)); it->second = findIt->second; changed = true; } @@ -145,14 +148,16 @@ namespace AMDiS { while (elInfo) { for (int i = 0; i < mesh->getGeo(VERTEX); i++) if (mapDelDofs.count(elInfo->getElement()->getDof(i)) == 1) - elInfo->getElement()->setDof(i, const_cast<DegreeOfFreedom*>(mapDelDofs[elInfo->getElement()->getDof(i)])); + elInfo->getElement()->setDof(i, const_cast<DegreeOfFreedom*>(mapDelDofs[elInfo->getElement()->getDof(i)])); elInfo = stack.traverseNext(elInfo); } + for (std::map<const DegreeOfFreedom*, const DegreeOfFreedom*>::iterator it = mapDelDofs.begin(); it != mapDelDofs.end(); ++it) mesh->freeDof(const_cast<DegreeOfFreedom*>(it->first), VERTEX); } + } diff --git a/AMDiS/src/parallel/MeshManipulation.h b/AMDiS/src/parallel/MeshManipulation.h index 03ec50b0b89cee8f5e36e8e6c6800c09b0a84a5d..2602f9591a6ef770ff8660a8006f6fffd20c2a7c 100644 --- a/AMDiS/src/parallel/MeshManipulation.h +++ b/AMDiS/src/parallel/MeshManipulation.h @@ -27,23 +27,22 @@ #include <map> #include "AMDiS_fwd.h" +#include "parallel/ElementObjectData.h" namespace AMDiS { class MeshManipulation { public: - MeshManipulation(Mesh *m) - : mesh(m) + MeshManipulation(FiniteElemSpace *space) + : feSpace(space), + mesh(space->getMesh()) {} - void deleteDoubleDofs(std::set<MacroElement*>& newMacroEl); - - private: - void deleteDoubleDofs(std::map<int, MacroElement*>& leafInMacroEl, - std::set<MacroElement*>& newMacroEl, - int mode); + void deleteDoubleDofs(std::set<MacroElement*>& newMacroEl, ElementObjects &elObj); private: + FiniteElemSpace *feSpace; + Mesh *mesh; }; diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc index f44144ceeeb1b724d835dac85f104c124b9d6907..f4d654268c831974db8f733e27b28ba18ac24222 100644 --- a/AMDiS/src/parallel/ParallelDebug.cc +++ b/AMDiS/src/parallel/ParallelDebug.cc @@ -288,7 +288,7 @@ namespace AMDiS { DOFIterator<WorldVector<double> > it(&coords, USED_DOFS); for (it.reset(); !it.end(); ++it) - coordsToIndex[(*it)] = pdb.mapLocalGlobalDofs[it.getDOFIndex()]; + coordsToIndex[(*it)] = pdb.mapLocalGlobalDofs[it.getDOFIndex()]; StdMpi<CoordsIndexMap> stdMpi(pdb.mpiComm, true); for (RankToDofContainer::iterator it = pdb.sendDofs.begin(); @@ -402,6 +402,37 @@ namespace AMDiS { } + void ParallelDebug::testDoubleDofs(Mesh *mesh) + { + FUNCNAME("ParallelDebug::testDoubleDofs()"); + + std::map<WorldVector<double>, DegreeOfFreedom> cMap; + int foundError = 0; + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); + while (elInfo) { + for (int i = 0; i < 3; i++) { + WorldVector<double> &c = elInfo->getCoord(i); + if (cMap.count(c) == 0) { + cMap[c] = elInfo->getElement()->getDof(i, 0); + } else { + if (cMap[c] != elInfo->getElement()->getDof(i, 0)) { + MSG("[DBG] Found two DOFs %d and %d with the same coords %f %f!\n", + cMap[c], elInfo->getElement()->getDof(i, 0), c[0], c[1]); + foundError = 1; + } + } + } + + elInfo = stack.traverseNext(elInfo); + } + + mpi::globalAdd(foundError); + TEST_EXIT(foundError == 0)("Error found on at least one rank!\n"); + } + + void ParallelDebug::printMapLocalGlobal(MeshDistributor &pdb, int rank) { if (rank == -1 || pdb.mpiRank == rank) { diff --git a/AMDiS/src/parallel/ParallelDebug.h b/AMDiS/src/parallel/ParallelDebug.h index a9f472c0d5d16086ec2e0fd2638c85ff1aa6f448..3692119987e9ccb9501e6eb4292fedc135bb30a1 100644 --- a/AMDiS/src/parallel/ParallelDebug.h +++ b/AMDiS/src/parallel/ParallelDebug.h @@ -82,6 +82,9 @@ namespace AMDiS { RankToDofContainer &sendDofs, RankToDofContainer &recvDofs); + /// Tests if there are multiple DOFs in mesh with the same coords. + static void testDoubleDofs(Mesh *mesh); + /** \brief * This function is used for debugging only. It prints all information from * the local to global dof mapping, see \ref mapLocalGlobalDofs.