From 0c1e2d14b97642b6a03cdbbb40416e1cc41ae315 Mon Sep 17 00:00:00 2001 From: Thomas Witkowski <thomas.witkowski@gmx.de> Date: Fri, 19 Oct 2012 13:40:05 +0000 Subject: [PATCH] DOES NOT COMPILE, BUT I MUST MERGE TO WORK HARD ON WEEKEND.... --- AMDiS/src/parallel/MatrixNnzStructure.cc | 20 +- AMDiS/src/parallel/MeshDistributor.cc | 133 ++++---- AMDiS/src/parallel/MeshDistributor.h | 7 +- .../src/parallel/ParallelCoarseSpaceMatVec.h | 4 +- AMDiS/src/parallel/ParallelDebug.cc | 16 +- AMDiS/src/parallel/ParallelDofMapping.cc | 123 +++---- AMDiS/src/parallel/ParallelDofMapping.h | 300 +++++++++++++++--- AMDiS/src/parallel/PetscSolverFeti.h | 23 +- 8 files changed, 401 insertions(+), 225 deletions(-) diff --git a/AMDiS/src/parallel/MatrixNnzStructure.cc b/AMDiS/src/parallel/MatrixNnzStructure.cc index 1997e382..ba8d9cdd 100644 --- a/AMDiS/src/parallel/MatrixNnzStructure.cc +++ b/AMDiS/src/parallel/MatrixNnzStructure.cc @@ -122,8 +122,8 @@ namespace AMDiS { const FiniteElemSpace *rowFeSpace = dofMat->getRowFeSpace(); const FiniteElemSpace *colFeSpace = dofMat->getColFeSpace(); - if (rowDofMap.isDefinedFor(rowFeSpace) == false || - colDofMap.isDefinedFor(colFeSpace) == false) + if (rowDofMap.isDefinedFor(rowComp) == false || + colDofMap.isDefinedFor(colComp) == false) continue; // === Prepare MTL4 iterators for the current DOF matrix. === @@ -137,8 +137,8 @@ namespace AMDiS { // === Iterate on all DOFs of the row mapping. === - DofMap::iterator rowIt = rowDofMap[rowFeSpace].begin(); - DofMap::iterator rowEndIt = rowDofMap[rowFeSpace].end(); + DofMap::iterator rowIt = rowDofMap[rowComp].begin(); + DofMap::iterator rowEndIt = rowDofMap[colComp].end(); for (; rowIt != rowEndIt; ++rowIt) { // Go to the corresponding matrix row (note, both the mapping and the @@ -165,7 +165,7 @@ namespace AMDiS { perMap->fillAssociations(rowFeSpace, rowIt->second.global, elObjDb, perRowAsc); - if (localMatrix || rowDofMap[rowFeSpace].isRankDof(*cursor)) { + if (localMatrix || rowDofMap[rowComp].isRankDof(*cursor)) { // === The current row DOF is a rank DOF, so create the === // === corresponding nnz values directly on rank's nnz data. === @@ -178,7 +178,7 @@ namespace AMDiS { TEST_EXIT_DBG(localPetscRowIdx >= 0 && localPetscRowIdx < nRankRows) ("Should not happen! \n Debug info: DOF = %d globalRowIndx = %d petscRowIdx = %d localPetscRowIdx = %d rStart = %d compontens = %d from %d nRankRows = %d\n", *cursor, - rowDofMap[rowFeSpace][*cursor].global, + rowDofMap[rowComp][*cursor].global, petscRowIdx, localPetscRowIdx, rankStartRowIndex, @@ -190,7 +190,7 @@ namespace AMDiS { if (localMatrix) { for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) - if (colDofMap[colFeSpace].isSet(col(*icursor))) + if (colDofMap[colComp].isSet(col(*icursor))) dnnz[localPetscRowIdx]++; } else { MultiIndex colDofIndex; @@ -198,7 +198,7 @@ namespace AMDiS { // Traverse all non zero entries in this row. for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) { - if (colDofMap[colFeSpace].find(col(*icursor), colDofIndex) == false) + if (colDofMap[colComp].find(col(*icursor), colDofIndex) == false) continue; // Set of periodic row associations (is empty, if row DOF is not @@ -208,7 +208,7 @@ namespace AMDiS { perMap->fillAssociations(colFeSpace, colDofIndex.global, elObjDb, perColAsc); if (perColAsc.empty()) { - if (colDofMap[colFeSpace].isRankDof(col(*icursor))) + if (colDofMap[colComp].isRankDof(col(*icursor))) dnnz[localPetscRowIdx]++; else onnz[localPetscRowIdx]++; @@ -252,7 +252,7 @@ namespace AMDiS { for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor) { - if (colDofMap[colFeSpace].find(col(*icursor), colDofIndex) == false) + if (colDofMap[colComp].find(col(*icursor), colDofIndex) == false) continue; int petscColIdx = (colDofMap.isMatIndexFromGlobal() ? diff --git a/AMDiS/src/parallel/MeshDistributor.cc b/AMDiS/src/parallel/MeshDistributor.cc index e77fba20..c2974377 100644 --- a/AMDiS/src/parallel/MeshDistributor.cc +++ b/AMDiS/src/parallel/MeshDistributor.cc @@ -74,7 +74,7 @@ namespace AMDiS { : problemStat(0), initialized(false), name("parallel"), - feSpaces(0), + componentSpaces(0), mesh(NULL), refineManager(NULL), info(10), @@ -148,7 +148,7 @@ namespace AMDiS { TEST_EXIT(mpiSize > 1) ("Parallelization does not work with only one process!\n"); - TEST_EXIT(feSpaces.size() > 0) + TEST_EXIT(componentSpaces.size() > 0) ("No FE space has been defined for the mesh distributor!\n"); TEST_EXIT(mesh)("No mesh has been defined for the mesh distributor!\n"); @@ -158,18 +158,18 @@ namespace AMDiS { bool doNext = false; do { doNext = false; - for (unsigned int i = 0; i < feSpaces.size() - 1; i++) { - if (feSpaces[i]->getBasisFcts()->getDegree() > - feSpaces[i + 1]->getBasisFcts()->getDegree()) { - const FiniteElemSpace *tmp = feSpaces[i + 1]; - feSpaces[i + 1] = feSpaces[i]; - feSpaces[i] = tmp; + for (unsigned int i = 0; i < uniqueFeSpaces.size() - 1; i++) { + if (uniqueFeSpaces[i]->getBasisFcts()->getDegree() > + uniqueFeSpaces[i + 1]->getBasisFcts()->getDegree()) { + const FiniteElemSpace *tmp = uniqueFeSpaces[i + 1]; + uniqueFeSpaces[i + 1] = uniqueFeSpaces[i]; + uniqueFeSpaces[i] = tmp; doNext = true; } } } while (doNext); - elObjDb.setFeSpace(feSpaces[0]); + elObjDb.setFeSpace(uniqueFeSpaces[0]); // If required, create hierarchical mesh level structure. createMeshLevelStructure(); @@ -191,7 +191,8 @@ namespace AMDiS { elObjDb.setData(partitionMap, levelData); #if (DEBUG != 0) - ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], dofMap, + ParallelDebug::writeDebugFile(uniqueFeSpaces[uniqueFeSpaces.size() - 1], + dofMap, debugOutputDir + "mpi-dbg", "dat"); #endif @@ -315,7 +316,7 @@ namespace AMDiS { ParallelDebug::testInteriorBoundary(*this); ParallelDebug::followBoundary(*this); - debug::writeMesh(feSpaces[0], -1, debugOutputDir + "macro_mesh"); + debug::writeMesh(uniqueFeSpaces[0], -1, debugOutputDir + "macro_mesh"); MSG("Debug mode tests finished!\n"); #endif @@ -363,25 +364,12 @@ namespace AMDiS { TEST_EXIT_DBG(probStat->getFeSpaces().size()) ("No FE spaces in stationary problem!\n"); - // === Add FE spaces from stationary problem to mesh distributor. === + TEST_EXIT(componentSpaces.size() == 0) + ("Parallelization of coupled problems is deactived at the moment!\n"); - for (unsigned int i = 0; i < probStat->getFeSpaces().size(); i++) { - bool foundFeSpace = false; - for (unsigned int j = 0; j < feSpaces.size(); j++) { - if (feSpaces[j] == probStat->getFeSpaces()[i]) - foundFeSpace = true; - - TEST_EXIT(feSpaces[j]->getMesh() == probStat->getFeSpaces()[i]->getMesh()) - ("MeshDistributor does not yet support usage of different meshes!\n"); - } - - if (!foundFeSpace) - feSpaces.push_back(probStat->getFeSpaces()[i]); - } - - TEST_EXIT(feSpaces.size() > 0)("Should not happen!\n"); - - mesh = feSpaces[0]->getMesh(); + componentSpaces = probStat->getComponentFeSpaces(); + uniqueFeSpaces = probStat->getFeSpaces(); + mesh = uniqueFeSpaces[0]->getMesh(); info = probStat->getInfo(); switch (mesh->getDim()) { @@ -678,8 +666,8 @@ namespace AMDiS { DofEdge dofEdge1 = edge1.first->getEdge(edge1.second); WorldVector<double> c0, c1; - mesh->getDofIndexCoords(dofEdge0.first, feSpaces[0], c0); - mesh->getDofIndexCoords(dofEdge0.second, feSpaces[0], c1); + mesh->getDofIndexCoords(dofEdge0.first, uniqueFeSpaces[0], c0); + mesh->getDofIndexCoords(dofEdge0.second, uniqueFeSpaces[0], c1); MSG("FOUND EDGE %d/%d <-> %d/%d\n", edge0.first->getIndex(), edge0.second, @@ -976,7 +964,7 @@ namespace AMDiS { MPI::Wtime() - first); #if (DEBUG != 0) - debug::writeMesh(feSpaces[0], -1, debugOutputDir + "mesh"); + debug::writeMesh(uniqueFeSpaces[0], -1, debugOutputDir + "mesh"); #endif // Because the mesh has been changed, update the DOF numbering and mappings. @@ -1196,8 +1184,10 @@ namespace AMDiS { int rank = 0; SerUtil::deserialize(in, rank); - for (unsigned int j = 0; j < feSpaces.size(); j++) - deserialize(in, data[rank][feSpaces[j]], dofIndexMap[feSpaces[j]]); + for (unsigned int j = 0; j < componentSpaces.size(); j++) + deserialize(in, + data[rank][componentSpaces[j]], + dofIndexMap[componentSpaces[j]]); } } @@ -1274,7 +1264,8 @@ namespace AMDiS { Parameters::get("dbg->write part mesh", writePartMesh); if (writePartMesh > 0 && repartitioningCounter == 0) ParallelDebug::writePartitioningFile(debugOutputDir + "partitioning", - repartitioningCounter, feSpaces[0]); + repartitioningCounter, + uniqueFeSpaces[0]); repartitioningCounter++; @@ -1305,7 +1296,7 @@ namespace AMDiS { // === Run mesh partitioner to calculate a new mesh partitioning. === - partitioner->setLocalGlobalDofMap(&(dofMap[feSpaces[0]].getMap())); + partitioner->setLocalGlobalDofMap(&(dofMap[uniqueFeSpaces[0]].getMap())); bool partitioningSucceed = partitioner->partition(elemWeights, ADAPTIVE_REPART); if (!partitioningSucceed) { @@ -1482,11 +1473,11 @@ namespace AMDiS { // Note that also if there are no macros to be deleted, this function will // update the number of elements, vertices, etc. of the mesh. - mesh->removeMacroElements(deleteMacroElements, feSpaces); + mesh->removeMacroElements(deleteMacroElements, uniqueFeSpaces); // === Remove double DOFs. === MeshManipulation meshManipulation(mesh); - meshManipulation.deleteDoubleDofs(feSpaces, newMacroEl, elObjDb); + meshManipulation.deleteDoubleDofs(uniqueFeSpaces, newMacroEl, elObjDb); mesh->dofCompress(); partitioner->createPartitionMap(partitionMap); @@ -1521,7 +1512,8 @@ namespace AMDiS { MSG("AMDiS runs in debug mode, so make some test ...\n"); ParallelDebug::writePartitioningFile(debugOutputDir + "partitioning", - repartitioningCounter, feSpaces[0]); + repartitioningCounter, + uniqueFeSpaces[0]); ParallelDebug::testAllElements(*this); ParallelDebug::testDoubleDofs(mesh); ParallelDebug::testInteriorBoundary(*this); @@ -1567,11 +1559,11 @@ namespace AMDiS { // === Create DOF communicator. === - dofComm.init(0, levelData, feSpaces); + dofComm.init(0, levelData, uniqueFeSpaces); dofComm.create(intBoundary); if (levelData.getLevelNumber() > 1) { - dofCommSd.init(0, levelData, feSpaces); + dofCommSd.init(0, levelData, uniqueFeSpaces); dofCommSd.create(intBoundarySd); } @@ -1583,8 +1575,8 @@ namespace AMDiS { int nLevels = levelData.getLevelNumber(); boundaryDofInfo.resize(nLevels); - for (unsigned int i = 0; i < feSpaces.size(); i++) { - const FiniteElemSpace *feSpace = feSpaces[i]; + for (unsigned int i = 0; i < uniqueFeSpaces.size(); i++) { + const FiniteElemSpace *feSpace = uniqueFeSpaces[i]; for (int level = 0; level < nLevels; level++) { @@ -1636,7 +1628,7 @@ namespace AMDiS { if (partitioner->getElementInRank()[(*it)->getIndex()] == false) macrosToRemove.insert(*it); - mesh->removeMacroElements(macrosToRemove, feSpaces); + mesh->removeMacroElements(macrosToRemove, uniqueFeSpaces); } @@ -1657,13 +1649,13 @@ namespace AMDiS { int nLevels = levelData.getLevelNumber(); TEST_EXIT_DBG(nLevels >= 1)("Should not happen!\n"); - dofMap.init(levelData, feSpaces, feSpaces); + dofMap.init(levelData, componentSpaces, uniqueFeSpaces); dofMap.setMpiComm(levelData.getMpiComm(0), 0); dofMap.setDofComm(dofComm); dofMap.clear(); if (nLevels > 1) { - dofMapSd.init(levelData, feSpaces, feSpaces); + dofMapSd.init(levelData, componentSpaces, uniqueFeSpaces); dofMapSd.setMpiComm(levelData.getMpiComm(1), 1); dofMapSd.setDofComm(dofCommSd); dofMapSd.clear(); @@ -1671,13 +1663,13 @@ namespace AMDiS { createBoundaryDofs(); - for (unsigned int i = 0; i < feSpaces.size(); i++) - updateLocalGlobalNumbering(dofMap, dofComm, feSpaces[i]); + for (unsigned int i = 0; i < uniqueFeSpaces.size(); i++) + updateLocalGlobalNumbering(dofMap, dofComm, uniqueFeSpaces[i]); dofMap.update(); if (nLevels > 1) { - for (unsigned int i = 0; i < feSpaces.size(); i++) - updateLocalGlobalNumbering(dofMapSd, dofCommSd, feSpaces[i]); + for (unsigned int i = 0; i < uniqueFeSpaces.size(); i++) + updateLocalGlobalNumbering(dofMapSd, dofCommSd, uniqueFeSpaces[i]); dofMapSd.update(); } @@ -1691,27 +1683,28 @@ namespace AMDiS { MSG("------------- Debug information -------------\n"); MSG("| number of levels: %d\n", nLevels); - MSG("| number of FE spaces: %d\n", feSpaces.size()); + MSG("| number of FE spaces: %d\n", uniqueFeSpaces.size()); - for (unsigned int i = 0; i < feSpaces.size(); i++) { - MSG("| FE space = %d (pointer adr %p):\n", i, feSpaces[i]); - MSG("| nRankDofs = %d\n", dofMap[feSpaces[i]].nRankDofs); - MSG("| nOverallDofs = %d\n", dofMap[feSpaces[i]].nOverallDofs); - MSG("| rStartDofs = %d\n", dofMap[feSpaces[i]].rStartDofs); + for (unsigned int i = 0; i < uniqueFeSpaces.size(); i++) { + MSG("| FE space = %d (pointer adr %p):\n", i, uniqueFeSpaces[i]); + MSG("| nRankDofs = %d\n", dofMap[uniqueFeSpaces[i]].nRankDofs); + MSG("| nOverallDofs = %d\n", dofMap[uniqueFeSpaces[i]].nOverallDofs); + MSG("| rStartDofs = %d\n", dofMap[uniqueFeSpaces[i]].rStartDofs); } if (nLevels > 1) { - for (unsigned int i = 0; i < feSpaces.size(); i++) { + for (unsigned int i = 0; i < uniqueFeSpaces.size(); i++) { MSG("| FE space = %d:\n", i); - MSG("| nRankDofs = %d\n", dofMapSd[feSpaces[i]].nRankDofs); - MSG("| nOverallDofs = %d\n", dofMapSd[feSpaces[i]].nOverallDofs); - MSG("| rStartDofs = %d\n", dofMapSd[feSpaces[i]].rStartDofs); + MSG("| nRankDofs = %d\n", dofMapSd[uniqueFeSpaces[i]].nRankDofs); + MSG("| nOverallDofs = %d\n", dofMapSd[uniqueFeSpaces[i]].nOverallDofs); + MSG("| rStartDofs = %d\n", dofMapSd[uniqueFeSpaces[i]].rStartDofs); } } // debug::writeElementIndexMesh(mesh, debugOutputDir + "elementIndex-" + // lexical_cast<string>(mpiRank) + ".vtu"); - ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], dofMap, + ParallelDebug::writeDebugFile(uniqueFeSpaces[uniqueFeSpaces.size() - 1], + dofMap, debugOutputDir + "mpi-dbg", "dat"); debug::testSortedDofs(mesh, elMap); @@ -1722,14 +1715,16 @@ namespace AMDiS { ParallelDebug::testGlobalIndexByCoords(*this); } #else - for (unsigned int i = 0; i < feSpaces.size(); i++) + for (unsigned int i = 0; i < uniqueFeSpaces.size(); i++) MSG("FE space %d: nRankDofs = %d nOverallDofs = %d\n", - i, dofMap[feSpaces[i]].nRankDofs, dofMap[feSpaces[i]].nOverallDofs); + i, dofMap[uniqueFeSpaces[i]].nRankDofs, + dofMap[uniqueFeSpaces[i]].nOverallDofs); int tmp = 0; Parameters::get(name + "->write parallel debug file", tmp); if (tmp) - ParallelDebug::writeDebugFile(feSpaces[feSpaces.size() - 1], dofMap, + ParallelDebug::writeDebugFile(uniqueFeSpaces[uniqueFeSpaces.size() - 1], + dofMap, debugOutputDir + "mpi-dbg", "dat"); #endif @@ -1786,8 +1781,8 @@ namespace AMDiS { // MPI::COMM_WORLD.Barrier(); [TODO: CHANGE BECAUSE NOT ALL RANKS HAVE PERIODIC MAP!!!] double first = MPI::Wtime(); - for (unsigned int i = 0; i < feSpaces.size(); i++) - createPeriodicMap(feSpaces[i]); + for (unsigned int i = 0; i < uniqueFeSpaces.size(); i++) + createPeriodicMap(uniqueFeSpaces[i]); // MPI::COMM_WORLD.Barrier(); INFO(info, 8)("Creation of periodic mapping needed %.5f seconds\n", @@ -2024,14 +2019,14 @@ namespace AMDiS { // Create a map from DOF indices to the corresponding DOF pointers. map<const FiniteElemSpace*, map<int, const DegreeOfFreedom*> > dofIndexMap; - for (unsigned int i = 0; i < feSpaces.size(); i++) { - ElementDofIterator elDofIter(feSpaces[i]); + for (unsigned int i = 0; i < uniqueFeSpaces.size(); i++) { + ElementDofIterator elDofIter(uniqueFeSpaces[i]); TraverseStack stack; ElInfo *elInfo = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL); while (elInfo) { elDofIter.reset(elInfo->getElement()); do { - dofIndexMap[feSpaces[i]][elDofIter.getDof()] = elDofIter.getDofPtr(); + dofIndexMap[uniqueFeSpaces[i]][elDofIter.getDof()] = elDofIter.getDofPtr(); } while (elDofIter.next()); elInfo = stack.traverseNext(elInfo); diff --git a/AMDiS/src/parallel/MeshDistributor.h b/AMDiS/src/parallel/MeshDistributor.h index df2297cf..fcb4dc51 100644 --- a/AMDiS/src/parallel/MeshDistributor.h +++ b/AMDiS/src/parallel/MeshDistributor.h @@ -132,6 +132,7 @@ namespace AMDiS { return mesh; } +#if 0 /// Returns an FE space from \ref feSpaces. inline const FiniteElemSpace* getFeSpace(unsigned int i = 0) { @@ -147,6 +148,7 @@ namespace AMDiS { { return feSpaces; } +#endif /// Returns the DOF mapping object, \ref dofMap. inline ParallelDofMapping& getDofMap() @@ -486,7 +488,10 @@ namespace AMDiS { string name; /// Finite element spaces of the problem. - vector<const FiniteElemSpace*> feSpaces; + //vector<const FiniteElemSpace*> feSpaces; + vector<const FiniteElemSpace*> componentSpaces; + vector<const FiniteElemSpace*> uniqueFeSpaces; + /// Mesh of the problem. Mesh *mesh; diff --git a/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h b/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h index 843152cb..3ceb8e73 100644 --- a/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h +++ b/AMDiS/src/parallel/ParallelCoarseSpaceMatVec.h @@ -204,14 +204,12 @@ namespace AMDiS { * space is not really necessary. Rewrite this! * * \param[in] component Component number of the system. - * \param[in] feSpace Finite element space of the component. * \param[in] dof DOF index * * \return True, if the dof is a coarse space DOF in the component. * False otherwise. */ inline bool isCoarseSpace(int component, - const FiniteElemSpace *feSpace, DegreeOfFreedom dof) { FUNCNAME("ParallelCoarseSpaceMatVec::isCoarseSpace()"); @@ -222,7 +220,7 @@ namespace AMDiS { TEST_EXIT_DBG(coarseSpaceMap.count(component)) ("Component %d has no coarse space defined!\n", component); - return (*(coarseSpaceMap[component]))[feSpace].isSet(dof); + return (*(coarseSpaceMap[component]))[component].isSet(dof); } diff --git a/AMDiS/src/parallel/ParallelDebug.cc b/AMDiS/src/parallel/ParallelDebug.cc index cc51b18e..f62563ac 100644 --- a/AMDiS/src/parallel/ParallelDebug.cc +++ b/AMDiS/src/parallel/ParallelDebug.cc @@ -164,7 +164,7 @@ namespace AMDiS { perMap.periodicDofAssociations[feSpace].begin(); it != perMap.periodicDofAssociations[feSpace].end(); ++it) { WorldVector<double> c; - pdb.mesh->getDofIndexCoords(it->first, pdb.feSpaces[0], c); + pdb.mesh->getDofIndexCoords(it->first, pdb.uniqueFeSpaces[0], c); int nAssoc = it->second.size(); } @@ -329,7 +329,7 @@ namespace AMDiS { clock_t first = clock(); // Get FE space with basis functions of the highest degree - const FiniteElemSpace *feSpace = pdb.feSpaces[pdb.feSpaces.size() - 1]; + const FiniteElemSpace *feSpace = pdb.uniqueFeSpaces[pdb.uniqueFeSpaces.size() - 1]; int testCommonDofs = 1; Parameters::get("dbg->test common dofs", testCommonDofs); @@ -471,7 +471,7 @@ namespace AMDiS { FUNCNAME("ParallelDebug::testGlobalIndexByCoords()"); // Get FE space with basis functions of the highest degree - const FiniteElemSpace *feSpace = pdb.feSpaces[pdb.feSpaces.size() - 1]; + const FiniteElemSpace *feSpace = pdb.uniqueFeSpaces[pdb.uniqueFeSpaces.size() - 1]; DOFVector<WorldVector<double> > coords(feSpace, "tmp"); pdb.mesh->getDofIndexCoords(feSpace, coords); @@ -649,7 +649,7 @@ namespace AMDiS { #if 0 if (rank == -1 || pdb.mpiRank == rank) { - const FiniteElemSpace *feSpace = pdb.feSpaces[0]; + const FiniteElemSpace *feSpace = pdb.uniqueFeSpaces[0]; cout << "====== DOF MAP LOCAL -> GLOBAL ====== " << endl; @@ -700,7 +700,7 @@ namespace AMDiS { dofit != rankDofs.end(); ++dofit) { cout << " " << **dofit << endl; WorldVector<double> coords; - pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpaces[0], coords); + pdb.mesh->getDofIndexCoords(*dofit, pdb.uniqueFeSpaces[0], coords); coords.print(); } @@ -709,7 +709,7 @@ namespace AMDiS { dofit != rankAllDofs.end(); ++dofit) { cout << " " << **dofit << endl; WorldVector<double> coords; - pdb.mesh->getDofIndexCoords(*dofit, pdb.feSpaces[0], coords); + pdb.mesh->getDofIndexCoords(*dofit, pdb.uniqueFeSpaces[0], coords); coords.print(); } } @@ -886,13 +886,13 @@ namespace AMDiS { if (followThisBound(it->rankObj.elIndex, it->neighObj.elIndex)) debug::writeLocalElementDofs(pdb.mpiRank, it->rankObj.elIndex, - pdb.feSpaces[0]); + pdb.uniqueFeSpaces[0]); for (InteriorBoundary::iterator it(pdb.intBoundary.other); !it.end(); ++it) if (followThisBound(it->rankObj.elIndex, it->neighObj.elIndex)) debug::writeLocalElementDofs(pdb.mpiRank, it->rankObj.elIndex, - pdb.feSpaces[0]); + pdb.uniqueFeSpaces[0]); } diff --git a/AMDiS/src/parallel/ParallelDofMapping.cc b/AMDiS/src/parallel/ParallelDofMapping.cc index 99bb8826..d84f3378 100644 --- a/AMDiS/src/parallel/ParallelDofMapping.cc +++ b/AMDiS/src/parallel/ParallelDofMapping.cc @@ -38,7 +38,7 @@ namespace AMDiS { } - FeSpaceDofMap::FeSpaceDofMap(MeshLevelData* ld) + ComponentDofMap::ComponentDofMap(MeshLevelData* ld) : levelData(ld), dofComm(NULL), feSpace(NULL), @@ -49,7 +49,7 @@ namespace AMDiS { } - void FeSpaceDofMap::clear() + void ComponentDofMap::clear() { dofMap.clear(); @@ -62,9 +62,9 @@ namespace AMDiS { } - void FeSpaceDofMap::update() + void ComponentDofMap::update() { - FUNCNAME("FeSpaceDofMap::update()"); + FUNCNAME("ComponentDofMap::update()"); // === Compute local indices for all rank owned DOFs. === @@ -89,18 +89,18 @@ namespace AMDiS { } - void FeSpaceDofMap::computeGlobalMapping() + void ComponentDofMap::computeGlobalMapping() { - FUNCNAME("FeSpaceDofMap::computeGlobalMapping()"); + FUNCNAME("ComponentDofMap::computeGlobalMapping()"); for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) it->second.global = it->second.local + rStartDofs; } - void FeSpaceDofMap::computeNonLocalIndices() + void ComponentDofMap::computeNonLocalIndices() { - FUNCNAME("FeSpaceDofMap::computeNonLocalIndices()"); + FUNCNAME("ComponentDofMap::computeNonLocalIndices()"); TEST_EXIT_DBG(dofComm)("No DOF communicator defined!\n"); @@ -178,19 +178,9 @@ namespace AMDiS { FUNCNAME("ParallelDofMapping::init()"); levelData = &ldata; - feSpaces = fe; - feSpacesUnique = uniqueFe; isNonLocal = b; - - - // === Init the mapping for all different FE spaces. === - for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin(); - it != feSpacesUnique.end(); ++it) { - addFeSpace(*it); - data[*it].setNeedGlobalMapping(isNonLocal); - data[*it].setNonLocal(isNonLocal); - } + data->init(fe, uniqueFe, isNonLocal); } @@ -200,9 +190,8 @@ namespace AMDiS { TEST_EXIT_DBG(levelData)("No mesh level data object defined!\n"); - for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin(); - it != feSpacesUnique.end(); ++it) - data[*it].clear(); + for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next()) + data->clear(); nRankDofs = -1; nLocalDofs = -1; @@ -217,9 +206,8 @@ namespace AMDiS { mpiComm = m; meshLevel = l; - for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin(); - it != feSpacesUnique.end(); ++it) - data[*it].setMpiComm(m, l); + for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next()) + it->setMpiComm(m, l); } @@ -230,35 +218,19 @@ namespace AMDiS { dofComm = &dc; // Add the DOF communicator also to all FE space DOF mappings. - for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin(); - it != feSpacesUnique.end(); ++it) - data[*it].setDofComm(dc); + for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next()) + it->setDofComm(dc); } - void ParallelDofMapping::addFeSpace(const FiniteElemSpace* feSpace) - { - FUNCNAME("ParallelDofMapping::addFeSpace()"); - - if (data.count(feSpace)) - data.find(feSpace)->second.clear(); - else - data.insert(make_pair(feSpace, FeSpaceDofMap(levelData))); - - data.find(feSpace)->second.setFeSpace(feSpace); - } - - int ParallelDofMapping::computeRankDofs() { FUNCNAME("ParallelDofMapping::computeRankDofs()"); int result = 0; - for (unsigned int i = 0; i < feSpaces.size(); i++) { - TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n"); - result += data[feSpaces[i]].nRankDofs; - } - + for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next()) + result += it->nRankDofs; + return result; } @@ -268,11 +240,9 @@ namespace AMDiS { FUNCNAME("ParallelDofMapping::computeLocalDofs()"); int result = 0; - for (unsigned int i = 0; i < feSpaces.size(); i++) { - TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n"); - result += data[feSpaces[i]].nLocalDofs; - } - + for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next()) + result += it->nLocalDofs; + return result; } @@ -282,11 +252,9 @@ namespace AMDiS { FUNCNAME("ParallelDofMapping::computeOverallDofs()"); int result = 0; - for (unsigned int i = 0; i < feSpaces.size(); i++) { - TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n"); - result += data.find(feSpaces[i])->second.nOverallDofs; - } - + for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next()) + result += it->nOverallDofs; + return result; } @@ -296,11 +264,9 @@ namespace AMDiS { FUNCNAME("ParallelDofMapping::computeStartDofs()"); int result = 0; - for (unsigned int i = 0; i < feSpaces.size(); i++) { - TEST_EXIT_DBG(data.count(feSpaces[i]))("Should not happen!\n"); - result += data.find(feSpaces[i])->second.rStartDofs; - } - + for (ComponentIterator &it = data->getIteratorComponent(); !it.end(); it.next()) + result += it->rStartDofs; + return result; } @@ -310,9 +276,8 @@ namespace AMDiS { FUNCNAME("ParallelDofMapping::update()"); // First, update all FE space DOF mappings. - for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin(); - it != feSpacesUnique.end(); ++it) - data[*it].update(); + for (ComponentIterator &it = data->getIteratorData(); !it.end(); it.next()) + it->update(); // Compute all numbers from this mappings. nRankDofs = computeRankDofs(); @@ -329,13 +294,17 @@ namespace AMDiS { { FUNCNAME("ParallelDofMapping::update()"); - for (vector<const FiniteElemSpace*>::iterator it = fe.begin(); - it != fe.end(); ++it) - if (find(feSpacesUnique.begin(), feSpacesUnique.end(), *it) == - feSpacesUnique.end()) - ERROR_EXIT("Wrong FE space!\n"); + ERROR_EXIT("DAS MUSS ICH MIR MAL UEBERLEGEN!\n"); + + // for (vector<const FiniteElemSpace*>::iterator it = fe.begin(); + // it != fe.end(); ++it) + // if (find(feSpacesUnique.begin(), feSpacesUnique.end(), *it) == + // feSpacesUnique.end()) + // ERROR_EXIT("Wrong FE space!\n"); + + // feSpaces = fe; + - feSpaces = fe; update(); } @@ -357,9 +326,9 @@ namespace AMDiS { // Traverse all DOFs of the FE space and create for all rank owned DOFs // a matrix index. - DofMap& dofMap = data[feSpaces[i]].getMap(); + DofMap& dofMap = (*data)[i].getMap(); for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) { - if (data[feSpaces[i]].isRankDof(it->first)) { + if ((*data)[i].isRankDof(it->first)) { int globalMatIndex = it->second.local + offset; if (globalIndex) dofToMatIndex.add(i, it->second.global, globalMatIndex); @@ -370,7 +339,7 @@ namespace AMDiS { // Increase the offset for the next FE space by the number of DOFs owned // by the rank in the current FE space. - offset += data[feSpaces[i]].nRankDofs; + offset += (*data)[i].nRankDofs; // If there are no non local DOFs, continue with the next FE space. if (!isNonLocal) @@ -500,10 +469,10 @@ namespace AMDiS { ("No data for FE space at address %p!\n", feSpaces[firstComponent]); int firstRankDof = -1; - FeSpaceDofMap &feMap = data.find(feSpaces[firstComponent])->second; - DofMap &dofMap = feMap.getMap(); + ComponentDofMap &compMap = (*data)[firstComponent]; + DofMap &dofMap = compMap.getMap(); for (DofMap::iterator it = dofMap.begin(); it != dofMap.end(); ++it) { - if (feMap.isRankDof(it->first)) { + if (compMap.isRankDof(it->first)) { if (needMatIndexFromGlobal) firstRankDof = it->second.global; else @@ -517,7 +486,7 @@ namespace AMDiS { int nRankRows = 0; for (int i = firstComponent; i < firstComponent + nComponents; i++) - nRankRows += data.find(feSpaces[firstComponent])->second.nRankDofs; + nRankRows += (*data)[i].nRankDofs; ISCreateStride(mpiComm, nRankRows, firstMatIndex, 1, &is); } diff --git a/AMDiS/src/parallel/ParallelDofMapping.h b/AMDiS/src/parallel/ParallelDofMapping.h index da79892c..9a8e08d5 100644 --- a/AMDiS/src/parallel/ParallelDofMapping.h +++ b/AMDiS/src/parallel/ParallelDofMapping.h @@ -104,18 +104,18 @@ namespace AMDiS { * by the class \ref ParallelDofMapping to specifiy the mapping for a set of * FE spaces. */ - class FeSpaceDofMap + class ComponentDofMap { public: /// This constructor exists only to create std::map of this class and make /// use of the operator [] for read access. Should never be called. - FeSpaceDofMap() + ComponentDofMap() { ERROR_EXIT("Should not be called!\n"); } /// This is the only valid constructur to be used. - FeSpaceDofMap(MeshLevelData* ld); + ComponentDofMap(MeshLevelData* ld); /// Clears all data of the mapping. void clear(); @@ -152,7 +152,7 @@ namespace AMDiS { /// the rank. void insertRankDof(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1) { - FUNCNAME("FeSpaceDofMap::insertRankDof()"); + FUNCNAME("ComponentDofMap::insertRankDof()"); TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n"); @@ -166,7 +166,7 @@ namespace AMDiS { /// but is owned by a different rank, thus it is part of an interior boundary. void insertNonRankDof(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1) { - FUNCNAME("FeSpaceDofMap::insertNonRankDof()"); + FUNCNAME("ComponentDofMap::insertNonRankDof()"); TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n"); @@ -290,8 +290,234 @@ namespace AMDiS { /// int nRankDofs, nLocalDofs, nOverallDofs, rStartDofs; }; - - + + class ComponentIterator { + public: + virtual ComponentDofMap& operator*() = 0; + + virtual bool end() = 0; + + virtual void next() = 0; + + virtual void reset() = 0; + }; + + class ComponentDataInterface + { + public: + virtual ComponentDofMap& operator[](int compNumber) = 0; + + virtual ComponentDofMap& operator[](const FiniteElemSpace *feSpace) = 0; + + virtual bool isDefinedFor(int compNumber) const = 0; + + virtual ComponentIterator& getIteratorData() = 0; + + virtual ComponentIterator& getIteratorComponent() = 0; + + virtual void init(vector<const FiniteElemSpace*> &f0, + vector<const FiniteElemSpace*> &f1, + bool isNonLocal) = 0; + + protected: + /// The FE spaces for all components. + vector<const FiniteElemSpace*> feSpaces; + + /// The set of all FE spaces. It uniquly contains all different FE spaces + /// from \ref feSpaces. + vector<const FiniteElemSpace*> feSpacesUnique; + }; + + + class ComponentDataEqFeSpace : ComponentDataInterface + { + public: + ComponentDofMap& operator[](int compNumber) + { + const FiniteElemSpace *feSpace = feSpaces[compNumber]; + return componentData.find(feSpace)->second; + } + + ComponentDofMap& operator[](const FiniteElemSpace *feSpace) + { + TEST_EXIT_DBG(componentData.count(feSpace))("No data for FE space!\n");; + + return componentData.find(feSpace)->second; + } + + bool isDefinedFor(int compNumber) const + { + const FiniteElemSpace *feSpace = feSpaces[compNumber]; + return static_cast<bool>(componentData.count(feSpace)); + } + + ComponentIterator& getIteratorData() + { + iterData.reset(); + return iterData; + } + + ComponentIterator& getIteratorComponent() + { + iterComponent.reset(); + return iterComponent; + } + + void init(vector<const FiniteElemSpace*> &f0, + vector<const FiniteElemSpace*> &f1, + bool isNonLocal) + { + // === Init the mapping for all different FE spaces. === + + // for (vector<const FiniteElemSpace*>::iterator it = feSpacesUnique.begin(); + // it != feSpacesUnique.end(); ++it) { + // addFeSpace(*it); + // data[*it].setNeedGlobalMapping(isNonLocal); + // data[*it].setNonLocal(isNonLocal); + // } + + // void ParallelDofMapping::addFeSpace(const FiniteElemSpace* feSpace) + // { + // FUNCNAME("ParallelDofMapping::addFeSpace()"); + + // if (data.count(feSpace)) + // data.find(feSpace)->second.clear(); + // else + // data.insert(make_pair(feSpace, ComponentDofMap(levelData))); + + // data.find(feSpace)->second.setFeSpace(feSpace); + // } + + } + + + protected: + class IteratorData : public ComponentIterator { + public: + ComponentDofMap& operator*() + { + } + + bool end() + { + } + + void next() + { + } + + void reset() + { + } + }; + + class IteratorComponent : public ComponentIterator { + public: + ComponentDofMap& operator*() + { + } + + bool end() + { + } + + void next() + { + } + + void reset() + { + } + }; + + + map<const FiniteElemSpace*, ComponentDofMap> componentData; + + IteratorData iterData; + + IteratorComponent iterComponent; + }; + + + class ComponentDataDiffFeSpace : ComponentDataInterface + { + public: + ComponentDofMap& operator[](int compNumber) + { + TEST_EXIT_DBG(componentData.count(compNumber))("No data for component %d!\n", compNumber); + + return componentData.find(compNumber)->second; + } + + ComponentDofMap& operator[](const FiniteElemSpace *feSpace) + { + ERROR_EXIT("BLUB\n"); + } + + ComponentIterator& getIteratorData() + { + iterData.reset(); + return iterData; + } + + ComponentIterator& getIteratorComponent() + { + iterComponent.reset(); + return iterComponent; + } + + bool isDefinedFor(int compNumber) const + { + return (static_cast<unsigned int>(compNumber) < componentData.size()); + } + + protected: + class IteratorData : public ComponentIterator { + public: + ComponentDofMap& operator*() + { + } + + bool end() + { + } + + void next() + { + } + + void reset() + { + } + }; + + class IteratorComponent : public ComponentIterator { + public: + ComponentDofMap& operator*() + { + } + + bool end() + { + } + + void next() + { + } + + void reset() + { + } + }; + + map<unsigned int, ComponentDofMap> componentData; + + IteratorData iterData; + + IteratorComponent iterComponent; + }; + + /** * Implements the mapping from sets of distributed DOF indices to local and * global indices. The mapping works for a given set of FE spaces. Furthermore, @@ -305,7 +531,6 @@ namespace AMDiS { dofComm(NULL), isNonLocal(true), needMatIndexFromGlobal(false), - feSpaces(0), nRankDofs(1), nLocalDofs(1), nOverallDofs(1), @@ -374,28 +599,31 @@ namespace AMDiS { return needMatIndexFromGlobal; } - /// Access the DOF mapping for a given FE space. - inline FeSpaceDofMap& operator[](const FiniteElemSpace* feSpace) - { - FUNCNAME("ParallelDofMapping::operator[]()"); - - TEST_EXIT_DBG(data.count(feSpace)) - ("No data for FE space at address %p!\n", feSpace); - return data.find(feSpace)->second; + /// Access the DOF mapping for a given component number. + inline ComponentDofMap& operator[](int compNumber) + { + return (*data)[compNumber]; } - inline bool isDefinedFor(const FiniteElemSpace* feSpace) const + /// Access the DOF mapping for a given FE space + inline ComponentDofMap& operator[](const FiniteElemSpace *feSpace) { - return static_cast<bool>(data.count(feSpace)); + return (*data)[feSpace]; } - /// Returns the number of solution components the mapping is defined on. - inline int getNumberOfComponents() const + inline bool isDefinedFor(int compNumber) const { - return static_cast<int>(feSpaces.size()); + return data->isDefinedFor(compNumber); } + /// Returns the number of solution components the mapping is defined on. + // inline int getNumberOfComponents() const + // { + // ERROR_EXIT("WRITE SOMETHING MEANINGFUL!\n"); + // // return static_cast<int>(feSpaces.size()); + // } + /// Returns \ref nRankDofs, thus the number of DOFs owned by the rank. inline int getRankDofs() { @@ -453,19 +681,14 @@ namespace AMDiS { /// component number. inline int getLocalMatIndex(int ithComponent, DegreeOfFreedom d) { - FUNCNAME("ParallelDofMapping::getLocalMatIndex()"); - - TEST_EXIT_DBG(data[feSpaces[ithComponent]].isRankDof(d)) - ("Should not happen!\n"); - return dofToMatIndex.get(ithComponent, d) - rStartDofs; } - inline const FiniteElemSpace* getFeSpace(int i = 0) - { - TEST_EXIT_DBG(i < feSpacesUnique.size())("Wrong component number!\n"); - return feSpacesUnique[i]; - } + // inline const FiniteElemSpace* getFeSpace(int i = 0) + // { + // TEST_EXIT_DBG(i < feSpacesUnique.size())("Wrong component number!\n"); + // return feSpacesUnique[i]; + // } // Writes all data of this object to an output stream. void serialize(ostream &out) @@ -507,9 +730,6 @@ namespace AMDiS { } protected: - /// Insert a new FE space DOF mapping for a given FE space. - void addFeSpace(const FiniteElemSpace* feSpace); - /// Compute \ref nRankDofs. int computeRankDofs(); @@ -548,15 +768,8 @@ namespace AMDiS { /// mapping will be stored on global DOF indices. bool needMatIndexFromGlobal; - /// Maps from FE space pointers to DOF mappings. - map<const FiniteElemSpace*, FeSpaceDofMap> data; - - /// The FE spaces for all components. - vector<const FiniteElemSpace*> feSpaces; - - /// The set of all FE spaces. It uniquly contains all different FE spaces - /// from \ref feSpaces. - vector<const FiniteElemSpace*> feSpacesUnique; + /// Maps from components to DOF mappings. + ComponentDataInterface *data; /// Number of DOFs owned by rank. int nRankDofs; @@ -574,7 +787,6 @@ namespace AMDiS { /// consideration of possibly multiple components. DofToMatIndex dofToMatIndex; }; - } #endif diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h index c7035a86..6daa0cae 100644 --- a/AMDiS/src/parallel/PetscSolverFeti.h +++ b/AMDiS/src/parallel/PetscSolverFeti.h @@ -208,26 +208,23 @@ namespace AMDiS { void printStatistics(); - /// Checks whether a given DOF in a given FE space is a primal DOF. - inline bool isPrimal(const FiniteElemSpace *feSpace, - DegreeOfFreedom dof) + /// Checks whether a given DOF is a primal DOF in a given component. + inline bool isPrimal(int component, DegreeOfFreedom dof) { - return primalDofMap[feSpace].isSet(dof); + return primalDofMap[component].isSet(dof); } - /// Checks whether a given DOF in a give FE space is a dual DOF. - inline bool isDual(const FiniteElemSpace *feSpace, - DegreeOfFreedom dof) + /// Checks whether a given DOF is a dual DOF in a given component. + inline bool isDual(int component, DegreeOfFreedom dof) { - return dualDofMap[feSpace].isSet(dof); + return dualDofMap[component].isSet(dof); } - /// Checks whether a given DOF in a give FE space is an interface DOF. - inline bool isInterface(const FiniteElemSpace *feSpace, - DegreeOfFreedom dof) + /// Checks whether a given DOF is an interface DOF in a given component. + inline bool isInterface(int component, DegreeOfFreedom dof) { - if (feSpace == pressureFeSpace) - return interfaceDofMap[feSpace].isSet(dof); + if (component == pressureComponent) + return interfaceDofMap[component].isSet(dof); return false; } -- GitLab