diff --git a/AMDiS/src/CoarseningManager2d.cc b/AMDiS/src/CoarseningManager2d.cc index fb54607899b8202ed2a772b42c9229e7147f2ad9..0801118f41a22fdb3cf3da695ac67b5438cad251 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 c0087e7567368add236d5472d4e69ea9e2da9e26..a2405910248e6cb7cb807d92a65783986b264966 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 5993d2a5065a766a45237e87b54fc2b297f35903..de10fc950ae60a12cbb1dc929c7ed4f538e94977 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 27af395381c0a43e570b28b81f9f8943e3c90ad7..303bfcb13f1c9187decd41909cb3d237e8422afc 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 0ba49876a34f289dfad4fa7f9ad34a90e1c3cfa0..be0b9deb48cafe1f2e57f1554d30b65e2183b077 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 e4022eab75f1b13a41c0935aeb5b2ff3eeca8225..c47f82d4f945022105586d93f657fcd6b7afb6d5 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 eaa06e344c2a2c8eb05611b9a608737b75129c0a..59489b48ae1c4016292a9886a4aabddff22a2ad9 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 7617f52d6ec8ff1da9e908bc4245195f3ecda222..7e79326d04d6bd725db0f84706c5dfe2db36510a 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 2f243b29eaaf38d46caa6ca365e888a03be46794..8eef5d4388eecec1ce6a6c8599ef80ef19809971 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 1cd925470b4d902b58ac5978245ab3c8d424d54d..205e82be4f84e87073dd3550261c4c0892c24b3b 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 8f3d1cf967a1c003dad941c646c25e64cdb1fb87..e35e8dde6541422d38888723bf30f73592a84dd8 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 a0909ab6a7af4921b3baa49488923cfe10227b14..bb56836eae3875dd7e838ff6ba1f2df6f27b79fd 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;