diff --git a/AMDiS/src/Traverse.cc b/AMDiS/src/Traverse.cc index b8149f8db4afab787f00abff687b1aae94800900..ca1f20b1dae1642128d7daabad5811d4ed6dafc3 100644 --- a/AMDiS/src/Traverse.cc +++ b/AMDiS/src/Traverse.cc @@ -28,7 +28,7 @@ namespace AMDiS { ElInfo* TraverseStack::traverseFirst(Mesh *mesh, int level, Flag fill_flag) { FUNCNAME("TraverseStack::traverseFirst()"); - + TEST_EXIT_DBG(fill_flag.isSet(Mesh::CALL_REVERSE_MODE) == false || fill_flag.isSet(Mesh::CALL_EVERY_EL_PREORDER)) ("Not yet implemented!\n"); diff --git a/AMDiS/src/io/FileWriter.cc b/AMDiS/src/io/FileWriter.cc index ba198dc8e10a109579f8bb520b16f38630424c8c..b67440da15b1c1dba33f4579e37f43a834ced612 100644 --- a/AMDiS/src/io/FileWriter.cc +++ b/AMDiS/src/io/FileWriter.cc @@ -159,14 +159,6 @@ namespace AMDiS { // vtkVectorWriter.setWriteAs3dVector(writeAs3dVector); // vtkVectorWriter.writeFile(fn + paraviewFileExt); -#if HAVE_PARALLEL_DOMAIN_AMDIS - ERROR_EXIT("Simon must make some work in this, does not compile in parallel!\n"); -// if (MPI::COMM_WORLD.Get_rank() == 0) -// vtkVectorWriter.writeParallelFile(paraFilename + paraviewParallelFileExt, -// MPI::COMM_WORLD.Get_size(), -// filename, postfix); -#endif - MSG("ParaView file written to %s\n", (fn + paraviewFileExt).c_str()); } diff --git a/AMDiS/src/parallel/ElementObjectDatabase.cc b/AMDiS/src/parallel/ElementObjectDatabase.cc index a484363fa2d8a20f90baaba204ad3ec1f8a56fee..93308fb6937cc6621173c1f5f5507418076748d4 100644 --- a/AMDiS/src/parallel/ElementObjectDatabase.cc +++ b/AMDiS/src/parallel/ElementObjectDatabase.cc @@ -140,6 +140,12 @@ namespace AMDiS { if (periodicVertices.size() == 0) return; + // === Calculate smallest periodic boundary ID in mesh. === + + smallestPeriodicBcType = 0; + for (map<BoundaryType, VertexVector*>::iterator it = mesh->getPeriodicAssociations().begin(); + it != mesh->getPeriodicAssociations().end(); ++it) + smallestPeriodicBcType = std::min(smallestPeriodicBcType, it->first); // === Get all vertex DOFs that have multiple periodic associations. === diff --git a/AMDiS/src/parallel/ElementObjectDatabase.h b/AMDiS/src/parallel/ElementObjectDatabase.h index bb8d857ad7dfc067d9b71aa825b5df8765de3b7a..b68005ad799eac238b51ac7c7150dd239e7297d1 100644 --- a/AMDiS/src/parallel/ElementObjectDatabase.h +++ b/AMDiS/src/parallel/ElementObjectDatabase.h @@ -450,6 +450,14 @@ namespace AMDiS { return (periodicVertices.size() != 0); } + /// Returns true if the given boundary type is larger or equal to the smallest + /// periodic boundary ID in mesh. See \ref smallestPeriodicBcType for more + /// information. + bool isValidPeriodicType(BoundaryType t) const + { + return (t >= smallestPeriodicBcType); + } + /// Write the element database to disk. void serialize(ostream &out); @@ -574,15 +582,21 @@ namespace AMDiS { map<pair<BoundaryType, BoundaryType>, BoundaryType> bConnMap; - // The following three data structures store periodic DOFs, edges and faces. + /// The following three data structures store periodic DOFs, edges and faces. PerBoundMap<DegreeOfFreedom>::type periodicVertices; PerBoundMap<DofEdge>::type periodicEdges; PerBoundMap<DofFace>::type periodicFaces; - // Stores to each vertex all its periodic associations. + /// Defines the smallest boudary ID for periodic boundary conditions. This is + /// required to distinguish between "real" periodic boundaries and periodic + /// boundary IDs that are set by the parallel algorithm for indirectly + /// connected boundaries. + BoundaryType smallestPeriodicBcType; + + /// Stores to each vertex all its periodic associations. map<DegreeOfFreedom, std::set<BoundaryType> > periodicDofAssoc; - // Stores to each edge all its periodic associations. + /// Stores to each edge all its periodic associations. map<DofEdge, std::set<DofEdge> > periodicEdgeAssoc; map<pair<ElementObjectData, ElementObjectData>, bool> edgeReverseMode; diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index 452543424c4bb669a5c46c47734a7ea15513f035..c02706bef0a7881eabfd155750d0af5076f8c9fd 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -132,6 +132,8 @@ namespace AMDiS { { FUNCNAME("MeshDistributor::initParallelization()"); + MSG("RUN INIT\n"); + if (initialized) return; @@ -2067,13 +2069,6 @@ namespace AMDiS { { FUNCNAME("MeshDistributor::createPeriodicMap()"); - -// vector<ElementObjectData> objs = elObjDb.getElements(0); -// MSG("SIZE OF 0-DOF is: %d\n", objs.size()); - -// for (int i = 0; i < objs.size(); i++) -// MSG("0-DOF IN EL %d - %d\n", objs[i].elIndex, objs[i].ithObject); - StdMpi<vector<int> > stdMpi(mpiComm, false); // === Each rank traverse its periodic boundaries and sends the DOF === @@ -2112,15 +2107,8 @@ namespace AMDiS { DegreeOfFreedom globalDof1 = dofFeData[feSpace].mapDofToGlobal[*(dofs1[j])]; - if (!periodicMap.isPeriodicOnBound(feSpace, type, globalDof0)) { -// if (globalDof0 == 0) { -// MSG("A-ADD DOF 0 at BOUND: %d %d %d <-> %d %d %d \n", -// bound.rankObj.elIndex, bound.rankObj.subObj, bound.rankObj.ithObj, -// bound.neighObj.elIndex, bound.neighObj.subObj, bound.neighObj.ithObj); -// } - + if (!periodicMap.isPeriodicOnBound(feSpace, type, globalDof0)) periodicMap.add(feSpace, type, globalDof0, globalDof1); - } } } @@ -2132,9 +2120,6 @@ namespace AMDiS { for (vector<AtomicBoundary>::iterator boundIt = it->second.begin(); boundIt != it->second.end(); ++boundIt) { -// MSG(" NOW ADD DOFS ON PER BOUND EL %d %d with rank %d\n", -// boundIt->rankObj.elIndex, boundIt->rankObj.subObj, it->first); - int nDofs = dofs.size(); boundIt->rankObj.el->getAllDofs(feSpace, boundIt->rankObj, dofs); @@ -2177,14 +2162,8 @@ namespace AMDiS { // Check if this global DOF with the corresponding boundary type was // not added before by another periodic boundary from other rank. - if (!periodicMap.isPeriodicOnBound(feSpace, type, globalDofIndex)) { -// MSG("B-ADD DOF %d -> %d at BOUND from rank %d\n", -// globalDofIndex, -// mapGlobalDofIndex, -// it->first); - + if (!periodicMap.isPeriodicOnBound(feSpace, type, globalDofIndex)) periodicMap.add(feSpace, type, globalDofIndex, mapGlobalDofIndex); - } } } @@ -2201,7 +2180,9 @@ namespace AMDiS { for (vector<AtomicBoundary>::iterator boundIt = it->second.begin(); boundIt != it->second.end(); ++boundIt) { - if (boundIt->type >= -3) + // If the boundary type is a normal periodic BC, continue. We just + // operate on indirectly connected periodic boundaries. + if (elObjDb.isValidPeriodicType(boundIt->type)) continue; DofContainer dofs; @@ -2217,15 +2198,9 @@ namespace AMDiS { for (std::set<BoundaryType>::iterator perAscIt = assoc.begin(); perAscIt != assoc.end(); ++perAscIt) - if (*perAscIt >= -3) { -// MSG(" add map DOF %d -> %d (%d)\n", -// globalDof, -// periodicMap.map(feSpace, *perAscIt, globalDof), -// *perAscIt); - + if (elObjDb.isValidPeriodicType(*perAscIt)) perObjMap[*perAscIt][globalDof] = - periodicMap.map(feSpace, *perAscIt, globalDof); - } + periodicMap.map(feSpace, *perAscIt, globalDof); } } @@ -2240,7 +2215,7 @@ namespace AMDiS { // DOFs. I'm note sure what is correct here. for (map<int, PeriodicDofMap>::iterator it = stdMpi2.getRecvData().begin(); it != stdMpi2.getRecvData().end(); ++it) - periodicMap.add(feSpace, it->second); + periodicMap.add(feSpace, it->second, false); } diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h index 49c79d5da4a16453e3fe455a42a2fdbe0151a93d..95cea99cd2afcdaaff734e3dfb1f858f0c84aa6f 100644 --- a/AMDiS/src/parallel/MeshDistributor.h +++ b/AMDiS/src/parallel/MeshDistributor.h @@ -443,11 +443,13 @@ namespace AMDiS { void getAllBoundaryDofs(const FiniteElemSpace *feSpace, DofContainer& dofs); + const ElementObjectDatabase& getElementObjectDb() + { + return elObjDb; + } - public: /// Adds a stationary problem to the global mesh distributor objects. static void addProblemStatGlobal(ProblemStatSeq *probStat); - protected: void addProblemStat(ProblemStatSeq *probStat); diff --git a/AMDiS/src/parallel/PeriodicMap.cc b/AMDiS/src/parallel/PeriodicMap.cc index 31636ce06cf155fb0c919b65e482884e94c3d741..5dfd1b50873f032ee22384f422f4dad0e3bf1902 100644 --- a/AMDiS/src/parallel/PeriodicMap.cc +++ b/AMDiS/src/parallel/PeriodicMap.cc @@ -16,14 +16,15 @@ namespace AMDiS { void PeriodicMap::add(const FiniteElemSpace *feSpace, - PeriodicDofMap &newMap) + PeriodicDofMap &newMap, + bool a) { FUNCNAME("PeriodicMap::add()"); for (PeriodicDofMap::iterator it = newMap.begin(); it != newMap.end(); ++it) for (DofMap::iterator dofIt =it->second.begin(); dofIt != it->second.end(); ++dofIt) - add(feSpace, it->first, dofIt->second, dofIt->first); + add(feSpace, it->first, dofIt->second, dofIt->first, a); } diff --git a/AMDiS/src/parallel/PeriodicMap.h b/AMDiS/src/parallel/PeriodicMap.h index 11da11dab7c2ac0c5c93790870fc5761c88bc7a3..8a0814880f5a0bf59c7b449d202b90bd7375e325 100644 --- a/AMDiS/src/parallel/PeriodicMap.h +++ b/AMDiS/src/parallel/PeriodicMap.h @@ -101,7 +101,8 @@ namespace AMDiS { /// this DOFs that maps to a different DOF index than the given one. inline void add(const FiniteElemSpace *feSpace, BoundaryType type, - DegreeOfFreedom dof0, DegreeOfFreedom dof1) + DegreeOfFreedom dof0, DegreeOfFreedom dof1, + bool a = true) { FUNCNAME("PeriodicMap::map()"); @@ -110,12 +111,13 @@ namespace AMDiS { ("Should not happen!\n"); periodicDofMap[feSpace][type][dof0] = dof1; - periodicDofAssociations[feSpace][dof0].insert(type); + if (a) + periodicDofAssociations[feSpace][dof0].insert(type); } /// Adds a whole periodic mapping to the current one. - void add(const FiniteElemSpace *feSpace, PeriodicDofMap &newMap); + void add(const FiniteElemSpace *feSpace, PeriodicDofMap &newMap, bool a = true); /// For a given global DOF index, this function returns the set of periodic diff --git a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc index a621412e6a452635b1e86d9c0313ac510e546421..3e30747d17c7579d447c6e56175278d1a02a9ad8 100644 --- a/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc +++ b/AMDiS/src/parallel/PetscSolverGlobalMatrix.cc @@ -307,7 +307,7 @@ namespace AMDiS { std::set<int>& perColAsc = perMap.getAssociations(colFe, globalColDof); for (std::set<int>::iterator it = perColAsc.begin(); it != perColAsc.end(); ++it) - if (*it >= -3) + if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it)) perAsc.insert(*it); // Scale value to the number of periodic associations of the column index. @@ -322,7 +322,7 @@ namespace AMDiS { // First, add the original matrix index. newCols.push_back(globalColDof); - + // And add all periodic matrix indices. for (std::set<int>::iterator it = perAsc.begin(); it != perAsc.end(); ++it) { @@ -377,14 +377,14 @@ namespace AMDiS { std::set<int>& perColAsc = perMap.getAssociations(colFe, globalColDof); for (std::set<int>::iterator it = perColAsc.begin(); it != perColAsc.end(); ++it) - if (*it >= -3) + if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it)) perAsc.insert(*it); } std::set<int>& perRowAsc = perMap.getAssociations(rowFe, globalRowDof); for (std::set<int>::iterator it = perRowAsc.begin(); it != perRowAsc.end(); ++it) - if (*it >= -3) + if (meshDistributor->getElementObjectDb().isValidPeriodicType(*it)) perAsc.insert(*it); // Scale the value with respect to the number of periodic associations.