From 165d98152dbb7bf7506bc304c4e70a32e967b471 Mon Sep 17 00:00:00 2001 From: Thomas Witkowski <thomas.witkowski@gmx.de> Date: Fri, 11 Dec 2009 13:34:06 +0000 Subject: [PATCH] Several bugfixes in adaptivity within parallel code. --- AMDiS/src/CoarseningManager2d.cc | 72 ++++++++++++-------------------- AMDiS/src/DOFVector.cc | 7 +++- AMDiS/src/DOFVector.h | 10 +++-- AMDiS/src/DOFVector.hh | 5 +++ AMDiS/src/DirichletBC.cc | 2 +- AMDiS/src/Element.cc | 13 +++--- AMDiS/src/MeshStructure.cc | 22 ++++++---- AMDiS/src/MeshStructure.h | 4 +- AMDiS/src/ParallelDomainBase.cc | 56 ++++++++++++++++++++----- AMDiS/src/ParallelDomainVec.cc | 6 ++- AMDiS/src/ProblemVec.cc | 51 ++++++++++++---------- AMDiS/src/StdMpi.h | 12 ++++++ 12 files changed, 160 insertions(+), 100 deletions(-) diff --git a/AMDiS/src/CoarseningManager2d.cc b/AMDiS/src/CoarseningManager2d.cc index fb546078..0801118f 100644 --- a/AMDiS/src/CoarseningManager2d.cc +++ b/AMDiS/src/CoarseningManager2d.cc @@ -28,19 +28,14 @@ namespace AMDiS { TEST_EXIT_DBG(child[0]->getMark() < 0 && child[1]->getMark() < 0) ("element %d with children[%d,%d] must not be coarsend!\n", el->getIndex(), child[0]->getIndex(), child[1]->getIndex()); - - if (mesh->getNumberOfDOFs(EDGE)) { - /****************************************************************************/ - /* remove dof from common edge of child[0] and child[1] */ - /****************************************************************************/ + + // remove dof from common edge of child[0] and child[1] + if (mesh->getNumberOfDOFs(EDGE)) mesh->freeDOF(const_cast<int*>(child[0]->getDOF(4)), EDGE); - } + // remove dof from the barycenters of child[0] and child[1] if (mesh->getNumberOfDOFs(CENTER)) { - /****************************************************************************/ - /* remove dof from the barycenters of child[0] and child[1] */ - /****************************************************************************/ - int node = mesh->getNode(CENTER); + int node = mesh->getNode(CENTER); mesh->freeDOF(const_cast<int*>(child[0]->getDOF(node)), CENTER); mesh->freeDOF(const_cast<int*>(child[1]->getDOF(node)), CENTER); @@ -70,8 +65,10 @@ namespace AMDiS { int n_neigh, int bound) { - Triangle *el = dynamic_cast<Triangle*>(const_cast<Element*>(coarsenList->getElement(0))); - Triangle *neigh = dynamic_cast<Triangle*>(const_cast<Element*>(coarsenList->getElement(1))); + Triangle *el = + dynamic_cast<Triangle*>(const_cast<Element*>(coarsenList->getElement(0))); + Triangle *neigh = + dynamic_cast<Triangle*>(const_cast<Element*>(coarsenList->getElement(1))); DegreeOfFreedom *dof[3]; dof[0] = const_cast<int*>(el->getChild(0)->getDOF(2)); @@ -82,34 +79,28 @@ namespace AMDiS { dof[1] = dof[2] = 0; } - int node = mesh->getNode(EDGE); if (mesh->getNumberOfDOFs(EDGE)) { - /****************************************************************************/ - /* get new dof on el at the midpoint of the coarsening edge */ - /****************************************************************************/ + int node = mesh->getNode(EDGE); + // get new dof on el at the midpoint of the coarsening edge + if (!el->getDOF(node + 2)) { el->setDOF(node + 2, mesh->getDOF(EDGE)); - if (neigh) { - neigh->setDOF(node + 2, const_cast<int*>(el->getDOF(node+2))); - } + if (neigh) + neigh->setDOF(node + 2, const_cast<int*>(el->getDOF(node + 2))); } } - if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(CENTER)) { - coarsenList->addDOFParents(n_neigh); - } + if (mesh->getNumberOfDOFs(EDGE) || mesh->getNumberOfDOFs(CENTER)) + coarsenList->addDOFParents(n_neigh); + + // restrict dof vectors to the parents on the patch - /****************************************************************************/ - /* restrict dof vectors to the parents on the patch */ - /****************************************************************************/ int nrAdmin = mesh->getNumberOfDOFAdmin(); for (int iadmin = 0; iadmin < nrAdmin; iadmin++) { - std::list<DOFIndexedBase*>::iterator it; DOFAdmin* admin = const_cast<DOFAdmin*>(&mesh->getDOFAdmin(iadmin)); - std::list<DOFIndexedBase*>::iterator end = admin->endDOFIndexed(); - for (it = admin->beginDOFIndexed(); it != end; ++it) { - (*it)->coarseRestrict(*coarsenList, n_neigh); - } + for (std::list<DOFIndexedBase*>::iterator it = admin->beginDOFIndexed(); + it != admin->endDOFIndexed(); ++it) + (*it)->coarseRestrict(*coarsenList, n_neigh); } coarsenTriangle(el); @@ -117,9 +108,8 @@ namespace AMDiS { if (neigh) coarsenTriangle(neigh); - /****************************************************************************/ - /* now, remove those dofs in the coarcening edge */ - /****************************************************************************/ + // now, remove those dofs in the coarcening edge + mesh->freeDOF(dof[0], VERTEX); if (mesh->getNumberOfDOFs(EDGE)) { mesh->freeDOF(dof[1], EDGE); @@ -145,24 +135,18 @@ namespace AMDiS { return 0; // single leaves don't get coarsened if (el->getChild(0)->getMark() >= 0 || el->getChild(1)->getMark() >= 0) { - /****************************************************************************/ - /* one of the children must not be coarsend; return :-( */ - /****************************************************************************/ + // one of the children must not be coarsend; return :( el->setMark(0); return 0; } if (!el->getChild(0)->isLeaf() || !el->getChild(1)->isLeaf()) { - /****************************************************************************/ - /* one of the children is not a leaf element; try again later on */ - /****************************************************************************/ + // one of the children is not a leaf element; try again later on doMore = true; return 0; } - /****************************************************************************/ - /* give the refinement edge the right orientation */ - /****************************************************************************/ + // give the refinement edge the right orientation if (el->getDOF(0,0) < el->getDOF(1,0)) { edge[0] = const_cast<int*>(el->getDOF(0)); @@ -180,9 +164,7 @@ namespace AMDiS { coarse_list.setCoarsePatch(1, el_info->getOppVertex(2) == 2); } - /****************************************************************************/ - /* check wether we can coarsen the patch or not */ - /****************************************************************************/ + // check wether we can coarsen the patch or not // ========================================================================== // === check for periodic boundary ========================================== diff --git a/AMDiS/src/DOFVector.cc b/AMDiS/src/DOFVector.cc index c0087e75..a2405910 100644 --- a/AMDiS/src/DOFVector.cc +++ b/AMDiS/src/DOFVector.cc @@ -9,6 +9,9 @@ namespace AMDiS { template<> void DOFVector<double>::coarseRestrict(RCNeighbourList& list, int n) { + MSG("Coarsen operation \"%d\" in vector \"%s\"\n", + coarsenOperation, this->name.c_str()); + switch (coarsenOperation) { case NO_OPERATION: return; @@ -26,7 +29,9 @@ namespace AMDiS { (const_cast<BasisFunction*>(feSpace->getBasisFcts()))->coarseInter(this, &list, n); break; default: - ERROR_EXIT("invalid coarsen operation\n"); + std::cout << "ERROR WITH PTR = " << this << std::endl; + ERROR_EXIT("Invalid coarsen operation \"%d\" in vector \"%s\"\n", + coarsenOperation, this->name.c_str()); } } diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h index 5993d2a5..de10fc95 100644 --- a/AMDiS/src/DOFVector.h +++ b/AMDiS/src/DOFVector.h @@ -53,7 +53,11 @@ namespace AMDiS { elementVector(3), boundaryManager(NULL), nBasFcts(0) - {} + { +#ifdef HAVE_PARALLEL_DOMAIN_AMDIS + rankDofs = NULL; +#endif + } DOFVectorBase(const FiniteElemSpace *f, std::string n); @@ -213,12 +217,13 @@ namespace AMDiS { #ifdef HAVE_PARALLEL_DOMAIN_AMDIS inline void setRankDofs(std::map<DegreeOfFreedom, bool> &dofmap) { - // rankDofs = &dofmap; + rankDofs = &dofmap; } inline bool isRankDof(DegreeOfFreedom dof) { TEST_EXIT_DBG(rankDofs)("No rank dofs set!\n"); + return (*rankDofs)[dof]; } #endif @@ -267,7 +272,6 @@ namespace AMDiS { int dim; #ifdef HAVE_PARALLEL_DOMAIN_AMDIS - public: std::map<DegreeOfFreedom, bool> *rankDofs; #endif }; diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh index 27af3953..303bfcb1 100644 --- a/AMDiS/src/DOFVector.hh +++ b/AMDiS/src/DOFVector.hh @@ -630,6 +630,7 @@ namespace AMDiS { coarsenOperation = rhs.coarsenOperation; this->operators = rhs.operators; this->operatorFactor = rhs.operatorFactor; + if (rhs.boundaryManager) { if (this->boundaryManager) delete this->boundaryManager; @@ -639,6 +640,10 @@ namespace AMDiS { this->boundaryManager = NULL; } +#ifdef HAVE_PARALLEL_DOMAIN_AMDIS + this->rankDofs = rhs.rankDofs; +#endif + return *this; } diff --git a/AMDiS/src/DirichletBC.cc b/AMDiS/src/DirichletBC.cc index 0ba49876..be0b9deb 100644 --- a/AMDiS/src/DirichletBC.cc +++ b/AMDiS/src/DirichletBC.cc @@ -54,7 +54,7 @@ namespace AMDiS { for (int i = 0; i < nBasFcts; i++) { #ifdef HAVE_PARALLEL_DOMAIN_AMDIS - // if (vector->isRankDof(dofIndices[i])) + if (vector->isRankDof(dofIndices[i])) #endif if (localBound[i] == boundaryType) { if (f) { diff --git a/AMDiS/src/Element.cc b/AMDiS/src/Element.cc index e4022eab..c47f82d4 100644 --- a/AMDiS/src/Element.cc +++ b/AMDiS/src/Element.cc @@ -580,17 +580,16 @@ namespace AMDiS { int s1 = el->getSideOfChild(0, ithSide, elType); int s2 = el->getSideOfChild(1, ithSide, elType); - code.nextElement(); - - if (s1 != -1) + if (s1 != -1) { + code.nextElement(); fitElementToMeshCode(refineManager, code, el->getFirstChild(), s1, el->getChildType(elType)); - - code.nextElement(); - - if (s2 != -1) + } + if (s2 != -1) { + code.nextElement(); fitElementToMeshCode(refineManager, code, el->getSecondChild(), s2, el->getChildType(elType)); + } } } diff --git a/AMDiS/src/MeshStructure.cc b/AMDiS/src/MeshStructure.cc index eaa06e34..59489b48 100644 --- a/AMDiS/src/MeshStructure.cc +++ b/AMDiS/src/MeshStructure.cc @@ -53,18 +53,19 @@ namespace AMDiS { commit(); } - void MeshStructure::init(Element *el, int ithSide, int elType) + void MeshStructure::init(Element *el, int ithSide, int elType, bool reverseOrder) { FUNCNAME("MeshStructure::init()"); clear(); - addAlongSide(el, ithSide, elType); + addAlongSide(el, ithSide, elType, reverseOrder); commit(); } - void MeshStructure::addAlongSide(Element *el, int ithSide, int elType) + void MeshStructure::addAlongSide(Element *el, int ithSide, int elType, + bool reverseOrder) { insertElement(el->isLeaf()); @@ -72,10 +73,17 @@ namespace AMDiS { int s1 = el->getSideOfChild(0, ithSide, elType); int s2 = el->getSideOfChild(1, ithSide, elType); - if (s1 != -1) - addAlongSide(el->getFirstChild(), s1, el->getChildType(elType)); - if (s2 != -1) - addAlongSide(el->getSecondChild(), s2, el->getChildType(elType)); + if (!reverseOrder) { + if (s1 != -1) + addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), reverseOrder); + if (s2 != -1) + addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), reverseOrder); + } else { + if (s2 != -1) + addAlongSide(el->getSecondChild(), s2, el->getChildType(elType), reverseOrder); + if (s1 != -1) + addAlongSide(el->getFirstChild(), s1, el->getChildType(elType), reverseOrder); + } } } diff --git a/AMDiS/src/MeshStructure.h b/AMDiS/src/MeshStructure.h index 7617f52d..7e79326d 100644 --- a/AMDiS/src/MeshStructure.h +++ b/AMDiS/src/MeshStructure.h @@ -44,7 +44,7 @@ namespace AMDiS { /// Creates a mesh structure code from a mesh object by traversing it in preorder. void init(Mesh *mesh); - void init(Element *el, int ithSide, int elType); + void init(Element *el, int ithSide, int elType, bool reverseOrder); void init(const std::vector<unsigned long int>& initCode, int n) { @@ -137,7 +137,7 @@ namespace AMDiS { /// Insert a new element to the structure code. Is used by the init function. void insertElement(bool isLeaf); - void addAlongSide(Element *el, int ithSide, int elType); + void addAlongSide(Element *el, int ithSide, int elType, bool reverseOrder); /// Merges two mesh structure codes to one structure code. void merge(MeshStructure *structure1, diff --git a/AMDiS/src/ParallelDomainBase.cc b/AMDiS/src/ParallelDomainBase.cc index 2f243b29..8eef5d43 100644 --- a/AMDiS/src/ParallelDomainBase.cc +++ b/AMDiS/src/ParallelDomainBase.cc @@ -741,6 +741,7 @@ namespace AMDiS { if (mesh->getChangeIndex() == lastMeshChangeIndex) return; + MSG("MESH CHANGED!\n"); // === Create mesh structure codes for all ranks boundary elements. === @@ -755,7 +756,7 @@ namespace AMDiS { MeshStructure elCode; elCode.init(boundIt->rankObj.el, boundIt->rankObj.ithObj, - boundIt->rankObj.elType); + boundIt->rankObj.elType, false); sendCodes[it->first].push_back(elCode); } } @@ -763,11 +764,7 @@ namespace AMDiS { StdMpi<MeshCodeVec, MeshCodeVec> stdMpi(mpiComm); stdMpi.prepareCommunication(true); stdMpi.send(sendCodes); - - for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin(); - it != otherIntBoundary.boundary.end(); ++it) - stdMpi.recv(it->first); - + stdMpi.recv(otherIntBoundary.boundary); stdMpi.startCommunication<unsigned long int>(); @@ -775,6 +772,19 @@ namespace AMDiS { bool meshFitTogether = true; + if (mpiRank == 0) { + DOFVector<double> ddd(feSpace, "tmp"); + ddd.set(0.0); + VtkWriter::writeFile(&ddd, "testold0.vtu"); + } + + if (mpiRank == 1) { + DOFVector<double> ddd(feSpace, "tmp"); + ddd.set(0.0); + VtkWriter::writeFile(&ddd, "testold1.vtu"); + } + + for (RankToBoundMap::iterator it = otherIntBoundary.boundary.begin(); it != otherIntBoundary.boundary.end(); ++it) { @@ -784,15 +794,21 @@ namespace AMDiS { for (std::vector<AtomicBoundary>::iterator boundIt = it->second.begin(); boundIt != it->second.end(); ++boundIt) { + std::cout << "[" << mpiRank << "] GET CODE: "; + recvCodes[i].print(); + MeshStructure elCode; elCode.init(boundIt->rankObj.el, boundIt->rankObj.ithObj, - boundIt->rankObj.elType); + boundIt->rankObj.elType, true); + + std::cout << "[" << mpiRank << "] CALC CODE: "; + elCode.print(); if (elCode.getCode() != recvCodes[i].getCode()) { TEST_EXIT_DBG(refineManager)("Refinement manager is not set correctly!\n"); - // recvCodes[i].reset(); -// fitElementToMeshCode(refineManager, recvCodes[i], boundIt->rankObj.el, -// boundIt->rankObj.ithObj, boundIt->rankObj.elType); + recvCodes[i].reset(); + fitElementToMeshCode(refineManager, recvCodes[i], boundIt->rankObj.el, + boundIt->rankObj.ithObj, boundIt->rankObj.elType); meshFitTogether = false; } @@ -801,9 +817,27 @@ namespace AMDiS { } } + MSG("MESH FIT ALGO FINISHED!\n"); + + if (mpiRank == 0) { + DOFVector<double> ddd(feSpace, "tmp"); + ddd.set(0.0); + VtkWriter::writeFile(&ddd, "test0.vtu"); + } + + if (mpiRank == 1) { + DOFVector<double> ddd(feSpace, "tmp"); + ddd.set(0.0); + VtkWriter::writeFile(&ddd, "test1.vtu"); + } + + if (!meshFitTogether) { + MSG("MESH STRUCTURE CHANGED!\n"); + std::cout << "MESH HAS BEEN CHANGED!" << std::endl; - exit(0); + + exit(0); } updateLocalGlobalNumbering(); diff --git a/AMDiS/src/ParallelDomainVec.cc b/AMDiS/src/ParallelDomainVec.cc index 1cd92547..205e82be 100644 --- a/AMDiS/src/ParallelDomainVec.cc +++ b/AMDiS/src/ParallelDomainVec.cc @@ -83,6 +83,10 @@ namespace AMDiS { if (probVec->getRHS()->getDOFVector(i)->getBoundaryManager()) removeBoundaryCondition(const_cast<BoundaryIndexMap&>(probVec->getRHS()->getDOFVector(i)->getBoundaryManager()->getBoundaryConditionMap())); + + std::cout << "TEST SOL " << probVec->getSolution()->getDOFVector(i) << ": " << probVec->getSolution()->getDOFVector(i)->getCoarsenOperation() << std::endl; + + std::cout << "TEST RHS " << probVec->getRHS()->getDOFVector(i) << ": " << probVec->getRHS()->getDOFVector(i)->getCoarsenOperation() << std::endl; } // Remove periodic boundaries on elements in mesh. @@ -91,7 +95,7 @@ namespace AMDiS { while (elInfo) { elInfo->getElement()->deleteElementData(PERIODIC); elInfo = stack.traverseNext(elInfo); - } + } } void ParallelDomainVec::exitParallelization(AdaptInfo *adaptInfo) diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc index 8f3d1cf9..e35e8dde 100644 --- a/AMDiS/src/ProblemVec.cc +++ b/AMDiS/src/ProblemVec.cc @@ -338,7 +338,8 @@ namespace AMDiS { std::string numberedName = "rhs[" + boost::lexical_cast<std::string>(i) + "]"; rhs->setDOFVector(i, new DOFVector<double>(componentSpaces[i], numberedName)); - numberedName = name + boost::lexical_cast<std::string>(i); + + numberedName = "solution[" + boost::lexical_cast<std::string>(i) + "]"; solution->setDOFVector(i, new DOFVector<double>(componentSpaces[i], numberedName)); solution->getDOFVector(i)->setCoarsenOperation(COARSE_INTERPOL); @@ -546,7 +547,6 @@ namespace AMDiS { TIME_USED(first, clock())); #endif - } Flag ProblemVec::markElements(AdaptInfo *adaptInfo) @@ -558,13 +558,11 @@ namespace AMDiS { allowFirstRefinement(); Flag markFlag = 0; - for (int i = 0; i < nComponents; i++) { - if (marker[i]) { + for (int i = 0; i < nComponents; i++) + if (marker[i]) markFlag |= marker[i]->markMesh(adaptInfo, componentMeshes[i]); - } else { - WARNING("no marker for component %d\n", i); - } - } + else + WARNING("no marker for component %d\n", i); return markFlag; } @@ -813,6 +811,8 @@ namespace AMDiS { { FUNCNAME("ProblemVec::writeFiles()"); + std::cout << "TEST IN WRITE: " << solution->getDOFVector(0) << " : " << solution->getDOFVector(0)->getCoarsenOperation() << std::endl; + clock_t first = clock(); #ifdef _OPENMP @@ -1228,7 +1228,7 @@ namespace AMDiS { Mesh *mesh, Flag assembleFlag) { - /* ================ Initialization of vectors ==================== */ + // === Initialization of vectors === if (rhs->getBoundaryManager()) rhs->getBoundaryManager()->initVector(rhs); @@ -1236,15 +1236,11 @@ namespace AMDiS { solution->getBoundaryManager()->initVector(solution); #ifdef _OPENMP + // === Parallel boundary assemblage === + TraverseParallelStack stack; -#else - TraverseStack stack; -#endif - /* ================= Parallel Boundary Assemblage ================= */ -#ifdef _OPENMP #pragma omp parallel -#endif { // Each thread assembles on its own dof-vectors. DOFVector<double> *tmpRhsVec = new DOFVector<double>(rhs->getFESpace(), "tmpRhs"); @@ -1252,7 +1248,6 @@ namespace AMDiS { tmpRhsVec->set(0.0); tmpSolVec->set(0.0); - // (Parallel) traverse of mesh. ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag); while (elInfo) { @@ -1265,13 +1260,10 @@ namespace AMDiS { elInfo = stack.traverseNext(elInfo); } - // After (parallel) mesh traverse, the result is applied to the final // vectors. This section is not allowed to be executed by more than one // thread at the same time. -#ifdef _OPENMP #pragma omp critical -#endif { DOFVector<double>::Iterator rhsIt(rhs, USED_DOFS); DOFVector<double>::Iterator solIt(solution, USED_DOFS); @@ -1285,13 +1277,28 @@ namespace AMDiS { } } // pragma omp critical - delete tmpRhsVec; delete tmpSolVec; } // pragma omp parallel - + +#else + // === Sequentiell boundary assemblage === + + TraverseStack stack; + ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag); + while (elInfo) { + if (rhs->getBoundaryManager()) + rhs->getBoundaryManager()-> fillBoundaryConditions(elInfo, rhs); + + if (solution->getBoundaryManager()) + solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution); + + elInfo = stack.traverseNext(elInfo); + } +#endif + - /* ======================= Finalize vectors ================== */ + // === Finalize vectors === if (rhs->getBoundaryManager()) rhs->getBoundaryManager()->exitVector(rhs); diff --git a/AMDiS/src/StdMpi.h b/AMDiS/src/StdMpi.h index a0909ab6..bb56836e 100644 --- a/AMDiS/src/StdMpi.h +++ b/AMDiS/src/StdMpi.h @@ -149,6 +149,18 @@ namespace AMDiS { recvDataSize[fromRank] = size; } + template<class T> + void recv(std::map<int, T> &fromRanks) + { + FUNCNAME("StdMpi::recv()"); + + TEST_EXIT_DBG(commPrepared)("Communication is not prepared!\n"); + + for (typename std::map<int, T>::iterator it = fromRanks.begin(); + it != fromRanks.end(); ++it) + recvDataSize[it->first] = -1; + } + std::map<int, RecvT>& getRecvData() { return recvData; -- GitLab