diff --git a/AMDiS/src/DOFAdmin.cc b/AMDiS/src/DOFAdmin.cc index fa859f95e9750800d9888216b13bc4a5e23aa068..8d62aaedac798f56fd5c52f209cb07fa15a62fee 100644 --- a/AMDiS/src/DOFAdmin.cc +++ b/AMDiS/src/DOFAdmin.cc @@ -11,6 +11,7 @@ #include "DOFVector.h" #include "DOFIterator.h" #include "Serializer.h" +#include "VertexVector.h" namespace AMDiS { @@ -24,6 +25,7 @@ namespace AMDiS { init(); } + DOFAdmin::DOFAdmin(Mesh* m, std::string aName) : name(aName), mesh(m), @@ -33,9 +35,11 @@ namespace AMDiS { init(); } + DOFAdmin::~DOFAdmin() {} + void DOFAdmin::init() { firstHole = 0; @@ -46,6 +50,7 @@ namespace AMDiS { dofFree.clear(); } + DOFAdmin& DOFAdmin::operator=(const DOFAdmin& src) { if (this != &src) { @@ -67,10 +72,7 @@ namespace AMDiS { return *this; } - /****************************************************************************/ - /* use a bit vector to indicate used/unused dofs */ - /* storage needed: one bit per dof */ - /****************************************************************************/ + bool DOFAdmin::operator==(const DOFAdmin& ad) const { @@ -90,9 +92,10 @@ namespace AMDiS { ERROR_EXIT("TODO\n"); } - void DOFAdmin::freeDOFIndex(int dof) + + void DOFAdmin::freeDofIndex(int dof) { - FUNCNAME("DOFAdmin::freeDOFIndex()"); + FUNCNAME("DOFAdmin::freeDofIndex()"); TEST_EXIT_DBG(usedCount > 0)("no dofs in use\n"); TEST_EXIT_DBG((dof >= 0) && (dof < size))("invalid dof index %d\n",dof); @@ -107,7 +110,7 @@ namespace AMDiS { std::list<DOFContainer*>::iterator dcend = dofContainerList.end(); for (dc = dofContainerList.begin(); dc != dcend; ++dc) - (*dc)->freeDOFIndex(dof); + (*dc)->freeDofIndex(dof); dofFree[dof] = true; @@ -118,6 +121,7 @@ namespace AMDiS { holeCount++; } + int DOFAdmin::getDOFIndex() { FUNCNAME("DOFAdmin::getDOFIndex()"); @@ -157,6 +161,7 @@ namespace AMDiS { return dof; } + void DOFAdmin::enlargeDOFLists(int minsize) { FUNCNAME("DOFAdmin::enlargeDOFLists()"); @@ -184,6 +189,7 @@ namespace AMDiS { (*di)->resize(newval); } + void DOFAdmin::addDOFIndexed(DOFIndexedBase* dofIndexed) { FUNCNAME("DOFAdmin::addDOFIndexed()"); @@ -201,6 +207,7 @@ namespace AMDiS { } } + void DOFAdmin::removeDOFIndexed(DOFIndexedBase* dofIndexed) { FUNCNAME("DOFAdmin::removeDOFIndexed()"); @@ -221,16 +228,20 @@ namespace AMDiS { } } - TEST_EXIT(removed)("DOFIndexed not in list\n"); + TEST_EXIT(removed)("DOFIndexed not in list!\n"); } + void DOFAdmin::addDOFContainer(DOFContainer* cont) { FUNCNAME("DOFAdmin::addDOFContainer()"); + TEST_EXIT_DBG(cont)("no container\n"); + dofContainerList.push_back(cont); } + void DOFAdmin::removeDOFContainer(DOFContainer* cont) { FUNCNAME("DOFAdmin::removeDOFContainer()"); @@ -243,7 +254,8 @@ namespace AMDiS { return; } } - ERROR("container not in list\n"); + + ERROR("Container not in list!\n"); } @@ -262,7 +274,7 @@ namespace AMDiS { // mark used dofs DOFIteratorBase it(this, USED_DOFS); for (it.reset(); !it.end(); ++it) - new_dof[it.getDOFIndex()] = 1; + new_dof[it.getDOFIndex()] = 1; int n = 0, last = 0; for (int i = 0; i < size; i++) { /* create a MONOTONE compress */ @@ -288,7 +300,7 @@ namespace AMDiS { // get index of first changed dof int first = last; - for (int i = 0; i<size; i++) { + for (int i = 0; i < size; i++) { if (new_dof[i] < i && new_dof[i] >= 0) { first = i; break; @@ -303,8 +315,8 @@ namespace AMDiS { std::list<DOFContainer*>::iterator dc; std::list<DOFContainer*>::iterator endc = dofContainerList.end(); - for (dc = dofContainerList.begin(); dc != endc; dc++) - (*dc)->compressDOFContainer(n, new_dof); + for (dc = dofContainerList.begin(); dc != endc; ++dc) + (*dc)->compressDofContainer(n, new_dof); } @@ -317,6 +329,7 @@ namespace AMDiS { nrDOF[i] = v; } + void DOFAdmin::setNumberOfPreDOFs(int i, int v) { FUNCNAME("DOFAdmin::setNumberOfPreDOFs()"); @@ -326,11 +339,13 @@ namespace AMDiS { nr0DOF[i] = v; } + int DOFAdmin::calcMemoryUsage() { return sizeof(DOFAdmin); } + void DOFAdmin::serialize(std::ostream &out) { // write name @@ -352,7 +367,8 @@ namespace AMDiS { nrDOF.serialize(out); nr0DOF.serialize(out); -} + } + void DOFAdmin::deserialize(std::istream &in) { diff --git a/AMDiS/src/DOFAdmin.h b/AMDiS/src/DOFAdmin.h index 8214f3bd76f8a068d6d7b76b7851a6081e307f91..842b80c22bec5ce8a184166a88cff50b93966c0f 100644 --- a/AMDiS/src/DOFAdmin.h +++ b/AMDiS/src/DOFAdmin.h @@ -95,6 +95,11 @@ namespace AMDiS { /// Removes the given DOFContainer object from DOFAdmin. void removeDOFContainer(DOFContainer* dofContainer); + std::list<DOFContainer*>& getDCL() + { + return dofContainerList; + } + /** \brief * Removes all holes of unused DOF indices by compressing the used range of * indices (it does not resize the vectors). While the global index of a DOF @@ -257,7 +262,7 @@ namespace AMDiS { int getDOFIndex(); /// Frees index dof. Used by Mesh::getDOF() - void freeDOFIndex(int dof); + void freeDofIndex(int dof); /// void serialize(std::ostream &out); diff --git a/AMDiS/src/DOFContainer.h b/AMDiS/src/DOFContainer.h index 65d66bf80e252e5213236afa3e9311b9de7b3a57..b4349af96fb7bc5015d7654bac14bc5c9c8d1248 100644 --- a/AMDiS/src/DOFContainer.h +++ b/AMDiS/src/DOFContainer.h @@ -44,22 +44,25 @@ namespace AMDiS { */ virtual DegreeOfFreedom& operator[](int i) = 0; - virtual void freeDOFIndex(DegreeOfFreedom dof) {} + virtual void freeDofIndex(DegreeOfFreedom dof) {} /** \brief * Used by DOFAdmin to actualize the DOF indices in this container after * DOF compression. */ - virtual void compressDOFContainer(int size, std::vector<DegreeOfFreedom> &newDOF) + virtual void compressDofContainer(int size, std::vector<DegreeOfFreedom> &newDOF) { + FUNCNAME("DOFContainer::compressDofContainer()"); + for (int i = 0; i < size; i++) { int j = newDOF[operator[](i)]; if (j >= 0) operator[](i) = j; else - ERROR_EXIT("invalid dof in dof container\n"); + ERROR_EXIT("Invalid DOF %d in DOF container! (%d %d)\n", j, i, operator[](i)); } } + }; } diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h index cb20fdee4a47dc53e38c3f7f39ba7060b3990a2a..1e6037ec7d789859a0cf66c29d7aa4ef6132d2e2 100644 --- a/AMDiS/src/DOFVector.h +++ b/AMDiS/src/DOFVector.h @@ -447,7 +447,7 @@ namespace AMDiS { { FUNCNAME("DOFVector<T>::operator[]"); TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size())) - ("Illegal vector index %d.\n",i); + ("Illegal vector index %d.\n", i); return vec[i]; } @@ -457,7 +457,7 @@ namespace AMDiS { FUNCNAME("DOFVector<T>::operator[]"); TEST_EXIT_DBG(i >= 0 && i < static_cast<int>(vec.size())) - ("Illegal vector index %d.\n",i); + ("Illegal vector index %d.\n", i); return vec[i]; } @@ -639,7 +639,7 @@ namespace AMDiS { /** \brief * Implements DOFContainer::operator[]() by calling - * DOFVEctor<DegreeOfFreedom>::operator[]() + * DOFVector<DegreeOfFreedom>::operator[]() */ DegreeOfFreedom& operator[](DegreeOfFreedom i) { diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc index ce387cfaea469a3bd197dd53ee2893398abefb19..e88757979a3b6cb1b3b60666528377dc503cc82f 100644 --- a/AMDiS/src/Debug.cc +++ b/AMDiS/src/Debug.cc @@ -191,7 +191,6 @@ namespace AMDiS { } return NULL; - } @@ -486,6 +485,32 @@ namespace AMDiS { } + void printElementRefinementSequenz(Mesh *mesh, Element *el) + { + FUNCNAME("printElementRefinementSequenz()"); + + int elIndex = el->getIndex(); + std::vector<int> refSeq; + Element *parent = getParentElement(mesh, el); + + while (parent) { + if (parent->getChild(0) == el) + refSeq.push_back(0); + else + refSeq.push_back(1); + + el = parent; + parent = getParentElement(mesh, el); + } + + std::stringstream oss; + for (int i = refSeq.size() - 1; i >= 0; i--) + oss << refSeq[i]; + + MSG("Ref-Seq for element %d: %s\n", elIndex, oss.str().c_str()); + } + + int getLocalNeighbourIndex(Mesh *mesh, int elIndex, int neighIndex) { FUNCNAME("debug::getLocalNeighbourIndex()"); diff --git a/AMDiS/src/Debug.h b/AMDiS/src/Debug.h index e22caac8472eed9c9f0278321d1fd1f4a1ccb8a6..f383789d9ea84e8b74dbc09ffa815affb1c419e7 100644 --- a/AMDiS/src/Debug.h +++ b/AMDiS/src/Debug.h @@ -137,6 +137,8 @@ namespace AMDiS { void printElementHierarchie(Mesh *mesh, int elIndex); + void printElementRefinementSequenz(Mesh *mesh, Element *el); + int getLocalNeighbourIndex(Mesh *mesh, int elIndex, int neighIndex); /** \brief diff --git a/AMDiS/src/Mesh.cc b/AMDiS/src/Mesh.cc index d14d44a2dbe15cb71b2ac4111d9a5222b7e53c23..a71e89d2be70be47c146d18103fbbb03c7d55ed3 100644 --- a/AMDiS/src/Mesh.cc +++ b/AMDiS/src/Mesh.cc @@ -590,7 +590,7 @@ namespace AMDiS { TEST_EXIT_DBG(n + n0 <= ndof)("n=%d, n0=%d too large: ndof=%d\n", n, n0, ndof); for (int j = 0; j < n; j++) - localAdmin->freeDOFIndex(dof[n0 + j]); + localAdmin->freeDofIndex(dof[n0 + j]); } delete [] dof; diff --git a/AMDiS/src/RefinementManager2d.cc b/AMDiS/src/RefinementManager2d.cc index bbae784992708c41eb33fc86d4db29a55439b968..16126dffbdd32a80898bdb40082a3c658de4102d 100644 --- a/AMDiS/src/RefinementManager2d.cc +++ b/AMDiS/src/RefinementManager2d.cc @@ -51,7 +51,7 @@ namespace AMDiS { // === Get the refinement patch. === getRefinePatch(&elInfo, edge, 0, refineList, &n_neigh); - + // === Check for periodic boundary === diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc index b8085a50fe617be126aa905e51ba88e6b6884f6e..1a2e8c25da34dc52c6494e5be44b42c89e4b7fbe 100644 --- a/AMDiS/src/Traverse.cc +++ b/AMDiS/src/Traverse.cc @@ -740,6 +740,7 @@ namespace AMDiS { /* first, goto to leaf level, if necessary... */ if (!(el->isLeaf()) && neighbour < 2) { + if (stack_used >= static_cast<int>(elinfo_stack.size()) - 1) enlargeTraverseStack(); @@ -780,6 +781,7 @@ namespace AMDiS { while (stack_used > 1) { stack_used--; int elIsIthChild = info_stack[stack_used]; + if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE) && elIsIthChild != 0) elIsIthChild = (elIsIthChild == 1 ? 2 : 1); @@ -790,8 +792,8 @@ namespace AMDiS { nb = coarse_nb[elIsIthChild][nb]; if (nb == -1) - break; - + break; + TEST_EXIT_DBG(nb >= 0)("invalid coarse_nb %d\n",nb); } @@ -818,7 +820,7 @@ namespace AMDiS { if (traverse_mel == NULL) return NULL; nb = i; - + stack_used = 1; elinfo_stack[stack_used]->fillMacroInfo(const_cast<MacroElement*>(traverse_mel)); info_stack[stack_used] = 0; @@ -834,15 +836,27 @@ namespace AMDiS { if (stack_used >= stack_size - 1) enlargeTraverseStack(); - - int i = 2 - info_stack[stack_used]; - info_stack[stack_used] = i + 1; - int fillIthChild = i; - if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) + + TEST_EXIT_DBG(info_stack[stack_used] == 1 || info_stack[stack_used] == 2) + ("Should not happen!\n"); + + int fillIthChild = -1; + if (info_stack[stack_used] == 1) { + info_stack[stack_used] = 2; + fillIthChild = 1; + nb = 0; + } else { + info_stack[stack_used] = 1; + fillIthChild = 0; + nb = 1; + } + + if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) { fillIthChild = 1 - fillIthChild; + nb = 1 - nb; + } elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]); stack_used++; - nb = 1 - i; } /****************************************************************************/ @@ -861,12 +875,16 @@ namespace AMDiS { if (stack_used >= stack_size - 1) enlargeTraverseStack(); + if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) { + int t = 2 - nb; + info_stack[stack_used] = (t == 2 ? 1 : 2); + } else { + info_stack[stack_used] = 2 - nb; + } + int fillIthChild = 1 - nb; - if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) - fillIthChild = 1 - fillIthChild; elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]); - info_stack[stack_used] = 2 - nb; stack_used++; nb = 2; } @@ -882,12 +900,16 @@ namespace AMDiS { info_stack[stack_used] = i + 1; int fillIthChild = i; - if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) + nb = i; + + if (traverse_fill_flag.isSet(Mesh::CALL_REVERSE_MODE)) { fillIthChild = 1 - fillIthChild; + nb = 1 - i; + } + elinfo_stack[stack_used + 1]->fillElInfo(fillIthChild, elinfo_stack[stack_used]); stack_used++; - nb = i; elinfo = elinfo_stack[stack_used]; el = dynamic_cast<Triangle*>(const_cast<Element*>(elinfo->getElement())); diff --git a/AMDiS/src/VertexVector.h b/AMDiS/src/VertexVector.h index f44a587e14ecdc394bd9a1e92fcf12d4e3bf4729..8de5950e95881d4d5a01a456b0f835ee98500ca2 100644 --- a/AMDiS/src/VertexVector.h +++ b/AMDiS/src/VertexVector.h @@ -36,9 +36,9 @@ namespace AMDiS { void set(DegreeOfFreedom val); - void compressDOFContainer(int size, std::vector<DegreeOfFreedom> &newDOF) + void compressDofContainer(int size, std::vector<DegreeOfFreedom> &newDOF) { - DOFContainer::compressDOFContainer(size, newDOF); + DOFContainer::compressDofContainer(size, newDOF); int totalSize = getAdmin()->getSize(); for (int i = size; i < totalSize; i++) vec[i] = i; diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index 2f356178a0b5e3bc61c01701ed5513cf5d3b9e70..86d1790b34829c82006477683ebba552401d976c 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -425,6 +425,8 @@ namespace AMDiS { void MeshDistributor::removePeriodicBoundaryConditions() { + FUNCNAME("MeshDistributor::removePeriodicBoundaryConditions()"); + // Remove periodic boundaries in boundary manager on matrices and vectors. for (unsigned int i = 0; i < probStat.size(); i++) { int nComponents = probStat[i]->getNumComponents(); @@ -451,6 +453,12 @@ namespace AMDiS { elInfo->getElement()->deleteElementData(PERIODIC); elInfo = stack.traverseNext(elInfo); } + + // Remove periodic vertex associations + for (std::map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin(); + it != mesh->getPeriodicAssociations().end(); ++it) + const_cast<DOFAdmin&>(mesh->getDOFAdmin(0)).removeDOFContainer(dynamic_cast<DOFContainer*>(it->second)); + mesh->getPeriodicAssociations().clear(); } @@ -572,7 +580,7 @@ namespace AMDiS { bool meshFitTogether = true; for (RankToBoundMap::iterator it = allBound.begin(); it != allBound.end(); ++it) { - + MeshCodeVec &recvCodes = stdMpi.getRecvData()[it->first]; int i = 0; @@ -590,7 +598,7 @@ namespace AMDiS { boundIt->rankObj.subObj, boundIt->rankObj.ithObj, boundIt->rankObj.elType, - boundIt->neighObj.reverseMode); + boundIt->rankObj.reverseMode); if (b) meshFitTogether = false; @@ -640,7 +648,8 @@ namespace AMDiS { TEST_EXIT_DBG(elInfo->getElement() == el)("This should not happen!\n"); - meshChanged = fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode); + meshChanged = + fitElementToMeshCode2(code, stack, subObj, ithObj, elType, reverseMode); return meshChanged; } @@ -999,9 +1008,11 @@ namespace AMDiS { std::map<int, Element*> elIndexMap; std::map<int, int> elIndexTypeMap; + std::map<std::pair<DegreeOfFreedom, DegreeOfFreedom>, BoundaryType> periodicDofs; + std::map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc; + std::map<std::pair<DofEdge, DofEdge>, BoundaryType> periodicEdges; std::map<std::pair<DofFace, DofFace>, BoundaryType> periodicFaces; - std::map<std::pair<DegreeOfFreedom, DegreeOfFreedom>, BoundaryType> periodicDofs; // === PHASE 1 === @@ -1052,6 +1063,7 @@ namespace AMDiS { TEST_EXIT_DBG(neighDof > -1)("Should not happen!\n"); periodicDofs[std::make_pair(dof, neighDof)] = boundaryType; + periodicDofAssoc[dof].insert(boundaryType); } } } @@ -1076,6 +1088,41 @@ namespace AMDiS { elObjects.createRankData(); + + // === Search for interectly connected vertices in periodic boundaries. === + + if (periodicDofs.size() > 0) { + TEST_EXIT(mesh->getDim() == 2)("Parallel boundaries not yet supported in 3D!\n"); + + std::vector<DegreeOfFreedom> twoPeriodicDof; + for (std::map<DegreeOfFreedom, std::set<BoundaryType> >::iterator it = periodicDofAssoc.begin(); + it != periodicDofAssoc.end(); ++it) { + TEST_EXIT_DBG(it->second.size() < 3)("Should not happen!\n"); + + if (it->second.size() == 2) + twoPeriodicDof.push_back(it->first); + } + + TEST_EXIT_DBG(twoPeriodicDof.size() == 0 || twoPeriodicDof.size() == 4) + ("Should not happen (%d)!\n", twoPeriodicDof.size()); + + if (twoPeriodicDof.size() == 4) { + for (int i = 0; i < 4; i++) { + for (int j = 0; j < 4; j++) { + if (i == j) + continue; + + std::pair<DegreeOfFreedom, DegreeOfFreedom> perDofs = + std::make_pair(twoPeriodicDof[i], twoPeriodicDof[j]); + + if (periodicDofs.count(perDofs) == 0) + periodicDofs[perDofs] = 3; + } + } + } + } + + // === PHASE 3 === for (int geoPos = 0; geoPos < mesh->getDim(); geoPos++) { @@ -1813,75 +1860,6 @@ namespace AMDiS { } } - if (dofFromRank.size() > 0) { - TEST_EXIT_DBG(mesh->getDim() == 2) - ("Periodic boundary corner problem must be generalized to 3d!\n"); - } - - MPI::Request request[min(static_cast<int>(periodicBoundary.boundary.size() * 2), 4)]; - int requestCounter = 0; - std::vector<int*> sendBuffers, recvBuffers; - - for (std::map<DegreeOfFreedom, std::set<int> >::iterator it = dofFromRank.begin(); - it != dofFromRank.end(); ++it) { - if (it->second.size() == 2) { - TEST_EXIT_DBG(periodicDofAssociations[it->first].size() == 2) - ("DOF %d has only %d periodic associations!\n", - it->first, periodicDofAssociations[it->first].size()); - - int type0 = *(periodicDofAssociations[it->first].begin()); - int type1 = *(++(periodicDofAssociations[it->first].begin())); - - int *sendbuf = new int[2]; - sendbuf[0] = periodicDof[type0][it->first]; - sendbuf[1] = periodicDof[type1][it->first]; - - request[requestCounter++] = - mpiComm.Isend(sendbuf, 2, MPI_INT, *(it->second.begin()), 0); - request[requestCounter++] = - mpiComm.Isend(sendbuf, 2, MPI_INT, *(++(it->second.begin())), 0); - - sendBuffers.push_back(sendbuf); - - int *recvbuf1 = new int[2]; - int *recvbuf2 = new int[2]; - - request[requestCounter++] = - mpiComm.Irecv(recvbuf1, 2, MPI_INT, *(it->second.begin()), 0); - request[requestCounter++] = - mpiComm.Irecv(recvbuf2, 2, MPI_INT, *(++(it->second.begin())), 0); - - recvBuffers.push_back(recvbuf1); - recvBuffers.push_back(recvbuf2); - } - } - - MPI::Request::Waitall(requestCounter, request); - - int i = 0; - for (std::map<DegreeOfFreedom, std::set<int> >::iterator it = dofFromRank.begin(); - it != dofFromRank.end(); ++it) { - if (it->second.size() == 2) { - for (int k = 0; k < 2; k++) - for (int j = 0; j < 2; j++) - if (recvBuffers[i + k][j] != it->first) { - int globalDofIndex = it->first; - int mapGlobalDofIndex = recvBuffers[i + k][j]; - int type = 3; - - periodicDof[type][globalDofIndex] = mapGlobalDofIndex; - periodicDofAssociations[globalDofIndex].insert(type); - } - - i++; - } - } - - for (unsigned int i = 0; i < sendBuffers.size(); i++) - delete [] sendBuffers[i]; - - for (unsigned int i = 0; i < recvBuffers.size(); i++) - delete [] recvBuffers[i]; #if (DEBUG != 0) for (std::map<int, std::set<BoundaryType> >::iterator it = periodicDofAssociations.begin();